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ABSTRACT 


The  development  of  smart  structures  technology  has  coincided  with  the  increased  use  of 
composite  materials  in  structural  design.  Composite  laminates  have  forms  of  damage  that  are  not 
found  in  other  materials,  specifically  delamination  and  transverse  matrix  cracking.  An  in-depth 
understanding  of  the  effects  of  damage  on  smart  composite  structures  is  necessary  for  predicting 
not  only  the  life  of  the  structure,  but  also  for  modeling  any  method  to  be  used  for  damage 
detection.  The  objective  of  this  research  was  to  develop  a  comprehensive  model  for  accurately 
and  efficiently  modeling  smart  composite  structures  including  the  effects  of  composite  damage. 

First,  a  new,  efficient  method  for  modeling  smart  structures  with  piezoelectric  devices 
was  developed.  The  coupled  model  simultaneously  solves  for  the  mechanical  and  electrical 
response  of  the  system  using  mechanical  displacements  and  electrical  displacements.  The 
developed  theory  utilizes  a  refined  higher  order  displacement  field  that  accurately  captures  the 
transverse  shear  deformation  in  moderately  thick  laminates. 

The  model  was  then  extended  to  include  internal  damage  in  the  form  of  delamination  and 
matrix  cracking.  When  delamination  is  present,  the  sublaminates  are  modeled  as  individual 
plates  and  continuity  is  enforced  at  the  interfaces.  Matrix  cracking  was  modeled  as  a  reduction  in 
laminate  stiffness  using  parameters  determined  using  finite  element  analysis  of  a  representative 
crack. 

Finally,  the  simultaneous  optimization  of  both  mechanical  and  electrical  parameters  in  an 
adaptive  structural  system  was  studied.  This  study  demonstrates  how  multidisciplinary 
optimization  techniques,  such  as  the  Kreisselmeier-Steinhauser  function,  can  be  utilized  to 
optimize  both  structural  and  electrical  aspects  of  an  adaptive  structural  system.  Optimization  of 
piezoelectric  actuator  placement  and  electrical  circuitry  was  performed  on  passive  electrical 
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damping  systems. 
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Results  show  that  the  developed  model  is  capable  of  accurately  modeling  both  the 
mechanical  and  electrical  response  of  adaptive  structures.  These  results  show  that  traditional 
uncoupled  piezoelectric  modeling  techniques  do  not  take  into  account  many  electrical  effects, 
resulting  in  significant  errors  in  the  sensor  response  predicted  for  transient  loads.  The  developed 
model  provides  the  ability  to  predict  the  effect  of  composite  damage  on  the  behavior  of  adaptive 
structures  and  to  model  potential  methods  to  be  used  for  damage  detection. 
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1.  Introduction 


So-called  “smart”  or  “adaptive”  structures  are  structural  systems  that  utilize  distributed 
sensors  and  actuators  in  order  to  monitor  or  improve  structural  response.  A  common  use  of  the 
smart  structures  concept  is  the  use  of  piezoelectric  sensors,  actuators  and  a  control  system  in 
order  to  control  vibration  in  the  structure.  This  concept  can  also  be  used  to  reduce  the  noise  that  is 
transmitted  through  a  structure.  In  another  application  of  adaptive  structures  technology ,  smart 
materials  are  being  used  to  induce  shape  changes  in  structures,  such  as  wings  or  helicopter  rotors. 
This  would  allow  designers  to  be  able  to  control  the  lift  distribution  over  a  surface  and  improve 
aerodynamic  efficiency. 

The  development  of  smart  structures  technology  has  coincided  with  the  increased  use  of 
composite  materials  in  structural  design.  Composite  structures  are  becoming  increasingly  popular 
in  both  aerospace  and  other  systems  due  to  the  benefit  of  reduced  weight  for  given  strength  and 
stiffness  requirements.  However  composite  laminates  have  forms  of  damage  that  are  not  found  in 
other  materials,  specifically  delamination  and  transverse  matrix  cracking.  Both  of  these  forms  of 
damage  occur  when  the  composite  is  subjected  to  fatigue,  over-loading,  or  low-velocity  impact. 
These  forms  of  damage  do  not  cause  catastrophic  failure  of  the  structure,  but  they  do  affect  the 
way  in  which  it  responds  to  applied  loads.  Methods  are  being  developed  which  use  smart 
materials  to  detect  and  localize  damage  within  a  composite  structure.  However,  limitations  in 
modeling  techniques  of  adaptive  structures  has  caused  reliance  on  empirical  data  for  the 
development  of  damage  detection  methods. 

Modeling  of  smart  structures  has  traditionally  been  based  on  simplified  models  that 
generally  focus  only  on  the  aspects  that  are  of  interest  to  a  single  engineering  discipline. 
Structural  models  generally  treat  electronics  as  “black  boxes”  and  electronics  designers  generally 
ignore  the  structure  all  together.  However,  smart  structures  rely  on  materials  that  cross  many 
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engineering  disciplines  and  accurate  modeling  is  not  possible  unless  it  includes  the 
interdisciplinary  aspects  of  these  materials.  In  addition,  the  response  of  a  structure  is  not  constant 
over  time.  Environmental  changes  and  structural  damage  affect  the  response  of  a  smart  structural 
system.  An  in-depth  understanding  of  the  effects  of  damage  on  smart  composite  structures  is 
necessary  for  predicting  not  only  the  life  of  the  structure,  but  also  for  modeling  any  method  to  be 
used  for  damage  detection.  If  smart  materials  are  to  be  increasingly  used,  then  significant 
improvements  in  modeling  must  be  made  so  that  all  of  these  issues  can  be  addressed.  More 
accurate  models  will  improve  the  efficiency  of  the  control  system,  potentially  reducing  its  weight 
and  power  consumption. 

1.1.  Piezoelectric  materials 

Piezoelectric  materials  are  one  of  the  most  commonly  used  smart  materials.  Piezoelectric 
materials  create  electricity  when  pressure  is  applied  to  them  and  deform  when  exposed  to  an 
electric  field.  These  two  behaviors  are  known  as  the  direct  effect  and  converse  effect, 
respectively,  and  each  has  its  own  governing  equation.  There  are  both  ceramic  and  polymer 
forms  of  piezoelectric  materials.  The  most  commonly  used  piezoceramic  is  lead  zirconate 
titanate  (PZT).  Piezoceramics  can  generate  moderate  forces  and  have  excellent  frequency 
response,  but  are  heavier  than  most  structural  materials  and  are  rather  brittle.  Piezoelectric  films, 
such  as  Polyvinylidene  Flouride  (PVDF),  are  lightweight  and  flexible,  but  can  only  induce  small 
forces  in  a  structure. 

The  traditional  use  of  piezoelectric  materials  in  structural  applications  is  to  take  thin 
sheets  of  material  and  coat  each  side  with  an  electrically  conductive  metal  electrode.  This 
essentially  creates  a  capacitor  with  the  piezoelectric  material  acting  as  a  dielectric.  When  voltage 
is  applied  across  the  electrodes  an  electric  field  is  created  within  the  piezoelectric  material, 
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inducing  deformation.  When  the  piezoelectric  material  is  subjected  to  mechanical  strain,  an 
electric  field  is  induced  within  the  material  and  an  electrical  signal  is  generated.  Thus  the  same 
device  can  be  used  as  an  actuator  or  a  sensor  depending  on  the  electronics  connected  to  it.  By 
either  bonding  the  device  to  the  surface  of  a  structure  or  embedding  it  within  the  structure  it  is 
then  possible  to  apply  localized  forces  or  measure  local  strain. 

Piezoelectric  actuators  are  often  used  to  control  the  dynamic  response  of  structures, 
however  the  dynamic  response  of  the  structure  is  modified  by  the  addition  of  these  actuators.  The 
PZT  or  other  piezoelectric  material  adds  additional  stiffness  to  the  structure,  but  this  stiffness  is 
dependent  upon  the  electrical  circuitry  attached  to  each  PZT.  This  can  be  easily  seen  in  the 
simple  comparison  between  actuators  that  are  open  circuited  and  ones  which  are  shorted  out. 
Under  structural  deflection  potential  energy  is  stored  in  the  PZT  in  the  form  of  mechanical  strain 
and  also  electrical  potential  similar  to  a  capacitor.  When  the  PZT  is  shorted  out  by  having  its 
electrodes  connected,  electrical  potential  can  no  longer  be  stored.  This  reduction  in  the  systems 
ability  to  store  energy  results  in  a  reduction  in  the  effective  stiffness  of  the  system.  The 
magnitude  of  this  reduction  is  dependent  on  the  piezoelectric  material’s  dielectric  permittivity,  or 
capacitance,  relative  to  its  piezoelectric  properties.  In  more  practical  examples  this  coupling  is 
encountered  in  the  interaction  between  the  structural  stiffness  and  the  electrical  potential  resulting 
from  a  combination  of  mechanical  and  electrical  loading.  When  piezoelectric  constants  are 
experimentally  determined  the  PZT  is  usually  unconstrained.  When  applied  to  a  structure  the 
PZT  is  no  longer  able  to  deform  to  its  unconstrained  shape  when  a  voltage  is  applied.  This 
change  in  the  strain  field  causes  a  reduction  in  the  electrical  potential  stored  in  the  PZT  and,  thus, 
a  smaller  actuation  force  than  predicted  using  an  uncoupled  theory  based  on  the  unconstrained 
piezoelectric  constants.  For  this  reason  piezoelectric  manufacturers  often  list  free  and  clamped, 
as  well  as  shorted  and  open  circuited,  material  properties  for  PZT.  However  a  PZT  mounted  to  a 
flexible  structure  is  neither  free  nor  clamped,  so  arbitrary  use  of  the  manufacturer  values  can  lead 
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to  analytic  errors  during  structural  modeling.  If  the  piezoelectric-mechanical  coupling  were  to  be 
taken  into  account  then  these  separate  parameters  would  not  be  required. 

The  modeling  of  piezoelectric  materials  has  been  well  understood  for  some  time 
(Tiersten,  1969),  but  only  recently  have  researcher  began  to  model  structures  with  piezoelectric 
devices  attached  to  them.  Beams  with  surface  bonded  actuators  were  the  first  to  be  examined. 
These  beams  were  modeled  using  a  uniform  strain  model  by  Crawley  and  de  Luis  (1987),  and 
then  with  Euler-Bemoulli  beam  theory  by  Crawley  and  Anderson  (1989).  Twisting  of  beams 
using  piezoelectric  patches  was  modeled  by  Park,  Walz  and  Chopra  (1996). 

A  model  for  anisotropic  plates  with  piezoelectric  actuators  covering  the  entire  surface 
was  developed  by  Crawley  and  Lazarus  (1991).  Their  model  used  a  classical  plate  theory  and 
Rayleigh-Ritz  analysis.  Classical  plate  theory  was  used  by  Wang  and  Rogers  (1991)  to  determine 
the  equivalent  force  and  moment  induced  by  finite-length  piezoelectric  actuators  on  composite 
laminates.  A  finite  element  formulation  for  plates  with  distributed  actuators  and  sensors  was  later 

developed  by  Detwiler,  Shen  and  Venkayya  (1995). 

One  thing  that  all  of  these  methods  have  in  common  is  that  they  all  use  an  uncoupled 
approach  to  solving  for  the  structural  response.  In  an  uncoupled  approach  only  one  of  the 
governing  equations  associated  with  either  the  direct  effect  or  converse  effect  is  solved.  In 
modeling  an  actuator,  this  is  made  possible  by  assuming  that  the  electric  field  in  the  piezoelectric 
material  is  a  known  constant  value.  This,  of  course,  leads  to  two  separate  models  for  the 
structure,  depending  on  whether  the  piezoelectric  patch  is  being  used  as  an  actuator  or  as  a  sensor. 
It  also  leads  to  inaccuracies  in  the  results,  since  the  exact  transfer  of  energy  between  electrical 

and  mechanical  is  assumed  and  not  solved  for. 

Although  methods  using  the  coupled  piezoelectric  equations  have  been  known  for  a  long 
time  (Tiersten,  1967,  1969),  they  have  not  been  used  extensively  in  modeling  smart  structures.  In 
coupled  approaches  the  mechanical  and  electrical  response  (both  direct  and  converse  effects)  are 
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simultaneously  solved  for.  This  method  was  first  used  in  the  modeling  of  beams  by  Hagood, 
Chung  and  Flotow  (1990).  A  general  two-way  coupled  theory  was  introduced  in  the  hybrid  plate 
theory  developed  by  Mitchell  and  Reddy  (1995).  In  this  model  a  mathematical  description  of  the 
electric  potential  was  assumed  and  could  be  solved  for  simultaneous  with  mechanical 
displacements.  This  model  was  extended  to  a  refined  higher-order  plate  theory  and  coupled  to  the 
thermal  field  by  Chattopadhyay,  Li  and  Gu  (1999).  The  description  of  the  electrical  potential  was 
then  improved  by  using  a  third  order  function  of  the  through-the-thickness  coordinate  Zhou, 
Chattopadhyay  and  Gu  (2000).  This  model  allowed  for  conservation  of  charge  across  the  PZT 
and  showed  the  interaction  of  bending  strains  with  the  electrical  field  within  the  piezoelectric 
material.  The  disadvantage  of  these  models  is  that  additional  degrees  of  freedom  are  added  to  the 
system  to  describe  the  distribution  of  electrical  potential  (and  thus  the  electric  field)  within  each 
piezoelectric  device.  The  number  of  degrees  of  freedom  increases  with  the  PZT  area  and  level  of 
finite  element  mesh  refinement.  Also,  the  resulting  system  of  equations  is  no  longer  symmetric 
positive  definite,  which  adds  to  the  computational  time. 

A  description  of  the  electrical  interaction  with  the  structure  also  allows  numerous  other 
applications  such  as  modeling  of  passive  electrical  damping  systems  and  self-sensing  actuators, 
as  well  as  analytical  determination  of  power  consumption.  Passive  electrical  damping  is  a 
technique  in  which  a  tuned  electrical  circuit  is  used  to  damp  vibration  in  a  structure  by  dissipating 
energy  as  heat  in  a  resistor.  Optimal  circuit  parameters  for  damping  of  a  particular  vibrational 
mode  were  determined  by  Hagood  and  Flotow  (1991).  Wu  (1996,  1998)  has  developed  an 
alternate  damping  circuit  and  addressed  damping  of  multiple  modes.  Kahn  and  Wang  (1994)  and 
Tsai  and  Wang  (1996,  1999)  have  demonstrated  combined  active-passive  control  circuits.  Others 
(Agnes,  1994;  Hagood  and  Crawley,  1991)  have  examined  the  application  of  passive  damping 
circuits  on  real  practical  structures.  Self-sensing  actuators  have  been  examined  by  Anderson  and 
Hagood  (1994)  using  coupled  equations  similar  to  those  used  in  this  work.  By  measuring  both 
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voltage  and  charge  flow  they  have  showed  that  it  is  possible  to  simultaneously  use  a  piezoelectric 
patch  as  an  actuator  and  a  sensor.  However,  all  of  these  elegant  electrical  models  relied  on  a 
comparatively  simple  description  of  the  strain  field  in  the  piezoelectric  material  and  analytical 
modeling  of  specific  structures  contained  noticeable  errors  as  a  result  of  the  approximations  used. 
If  accurate  structural  models  could  be  coupled  with  these  electrical  concepts  then  very  useful  and 
accurate  design  tools  could  be  developed. 

1.2.  Delamination 

Delamination  occurs  when  a  region  of  composite  material  becomes  debonded  from  the 
layers  below  it.  A  considerable  amount  of  research  effort  has  been  expended  in  the  understanding 
and  modeling  of  this  failure  process.  However,  the  vast  majority  of  this  work  has  focused  on  the 
initiation  and  propagation  of  delamination  under  static  and  fatigue  loading,  and  also  the  effect  on 
the  laminate’s  load  carrying  capability.  Delamination  also  affects  the  dynamic  characteristics  of  a 
composite  structure,  but  these  affects  are  subtler  and  have  been  less  studied. 

A  number  of  investigations  into  modeling  vibration  in  delaminated  composite  laminates 
have  been  made.  Delamination  of  a  composite  beam  was  first  modeled  by  Ramkumar,  Kulkarm 
and  Pipes  (1979),  but  their  model  predicted  a  severe  drop  in  the  natural  frequencies  which  was 
inconsistent  with  experimental  results.  An  experimental  investigation  on  the  effects  of 
delamination  on  vibration  properties  was  performed  by  Lee,  Sun  and  Liu  (1987)  and  Grady  and 
Meyn  (1989).  Beams  were  further  studied  using  classical  laminate  theory,  including  bending- 
extension  coupling,  by  a  number  of  other  researchers  (Wang,  et  al.,  1982;  Mujumdar  and 
Suryanarayan,  1988;  Tracy  and  Pardoen,  1989;  Shen  and  Grady,  1992).  Composite  plates  have 
also  been  modeled,  first  by  Yin  and  Jane  (1988)  and  later  by  a  number  of  other  authors  (Zak,  et 
al.,  2000; ).  Plates  with  multiple  delaminations  have  been  modeled  by  Ju,  et  al.  (1994)  and  Wang 
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and  Lin  (1995).  The  effect  of  local  thickening  was  studied  by  Hou  and  Jeronimidis  (1999).  Most 
recently,  nonlinear  vibration  was  studied  by  Lu,  Lestari  and  Hanagud  (2001). 

Most  of  these  methods  used  classical  plate  theory,  but  others  have  used  more  complex 
models.  Transverse  shear  has  been  included  by  a  number  of  authors  (Kardomateas  and 
Schmueser,  1988;  Pavier  and  Clarke,  1996).  Three-dimensional  models  (Yang  and  He,  1994) 
have  been  shown  to  provide  accurate  results,  but  with  computational  expense.  Layer-wise 
theories  (Barbero  and  Reddy,  1991)  have  also  been  used,  but  the  computational  effort  can  be 
quite  large  for  thick  laminates  with  a  large  number  of  layers.  A  refined  higher  order  theory  was 
developed  by  Chattopadhyay  and  Gu  (1994),  and  it  has  been  shown  to  be  efficient  and  provide 

accurate  solutions  (Chattopadhyay  and  Gu,  1997,  1999). 

The  effect  of  contact  between  the  sublaminates  during  plate  vibration  has  been  ignored  in 
most  works.  This  contact  will  alter  the  local  response  of  the  structure  in  the  region  of  the 
delamination,  particularly  at  high  frequencies,  as  shown  by  Kwon  and  Aygunes  (1996),  but  the 
effect  on  the  global  response  has  not  been  studied.  The  contact  between  the  sublaminates  makes 
the  response  a  nonlinear  problem.  This  nonlinearity  can  be  overcome  by  using  time-integration 
methods.  However,  the  sudden  impact  of  the  sublaminates  creates  oscillations  in  analytical 
response  that  will  destabilize  most  numerical  methods.  Kwon  and  Aygunes  used  the  method 
proposed  by  Hughes,  et  al.  (1976)  which  attempts  to  conserve  momentum  during  the  impact.  A 
discontinuous  time-integration  method  proposed  by  Cho  and  Kim  (1999a,b)  has  also  been  used  to 
treat  the  contact  between  the  sublaminates.  This  method  allows  discontinuity  to  be  incorporated 
into  time-integration  without  having  any  other  constraints  enforced  on  the  nodes  during  contact 
and  release,  and  it  also  has  been  shown  to  be  more  stable  than  other  methods. 

Delamination  modeling  in  adaptive  structures  has  received  surprisingly  little  attention. 
Authors  studying  damage  detection  have  generally  relied  on  experimental  results  rather  than 
accurate  laminate  models.  Cao,  et  al.  (1998)  investigated  the  effect  of  piezoelectric  actuators  on 
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delamination  growth.  Delamination  modeling  and  detection  using  the  refined  higher  order 
laminate  theory  has  also  been  studied  (Chattopadhyay,  Nam  and  Dragomir-Daescu,  1999, 
Chattopadhyay  and  Dragomir-Daescu,  1999).  There  exists  a  need  for  accurate  modeling 
techniques  for  delaminated  composites  that  includes  piezoelectric  modeling. 


1.3.  Matrix  cracking 

Matrix  cracking  is  a  common  form  of  damage  in  composite  laminates  and  results  from  a 
variety  of  conditions  including  low-velocity  impact,  fatigue,  and  excessive  loading.  Transverse 
matrix  cracking  occurs  in  a  composite  with  a  relatively  brittle  matrix,  such  as  epoxy.  The 
composite  fibers  fail  at  higher  strains  than  the  epoxy,  so  that  when  the  structure  is  over-loaded  or 
subjected  to  impact  cracks  initially  form  in  the  matrix.  These  cracks  generally  form  through  the 
width  of  each  lamina  along  the  fiber  direction.  Because  most  of  the  structural  load  is  carried  by 
the  fibers,  the  cracks  in  the  matrix  do  not  result  in  catastrophic  failure  of  the  structure.  However, 
the  matrix  cracks  do  cause  a  reduction  in  the  overall  stiffness  of  the  structure.  Matrix  cracks  can 
open  under  load,  thus  creating  an  increase  in  the  perceived  strains  of  a  laminate  and  a  reduction  in 
the  effective  stiffness.  In  engineering  design  it  is  often  assumed  that  once  matrix  cracking  has 
initiated,  the  lamina  with  cracks  is  no  longer  capable  of  carrying  load.  This  technique  is  known 
as  the  ply  discount  method.  However,  the  cracked  layer  is,  in  reality,  still  carrying  some  of  the 
load,  and  the  reduction  in  laminate  stiffness  will  be  a  function  of  crack  density.  Being  able  to 
accurately  describe  the  changes  in  structural  behavior  allows  designers  to  predict  the  effects  of 
localized  damage  on  static  and  dynamic  characteristics.  Also,  the  modeling  of  matrix  cracking  is 
necessary  for  the  development  of  smart  composite  structures  capable  of  detecting  damage. 

A  large  amount  of  research  has  been  conducted  on  developing  models  capable  of 
describing  the  effects  of  transverse  matrix  cracking  on  the  stiffness  of  composite  laminates. 
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Highsmith  and  his  co-authors  (1981,  1982)  were  among  the  first  to  examine  matrix  cracking. 
They  used  a  shear  lag  model  that  assumes  that  the  stresses  applied  the  laminate  transfer  to  the 
cracked  layer  from  adjacent  layers  via  shear  deformation  of  a  thin  boundary  layer,  call  the  shear 
transfer  region.  The  model  is  relatively  simple  and  gives  good  correlation  with  experimental  data 
provided  the  thickness  of  the  shear  transfer  region  can  be  determined.  Unfortunately  it  is  difficult 
to  systematically  determine  this  thickness  for  an  arbitrary  laminate.  This  model  was  later 
improved  by  Lim  and  Hong  (1989),  and  used  to  predict  shear  stiffiiess  by  Tsai  and  Daniel  (1992). 

Another  method  for  analysis  of  matrix  cracking  is  the  self-consistent  scheme  utilized  by 
Laws,  Dvorak  and  Hejazi  (1983)  and  Dvorak  (1985).  This  model  estimates  the  degraded  stiffness 
of  a  single  layer  with  stacked  cracks  based  on  the  reduction  for  a  body  with  randomly  distributed 
microcracks  or  inclusions  with  identical  shape  and  size.  The  resulting  degraded  stiffness 
components  for  each  layer  are  then  substituted  into  the  laminate  model  to  give  the  overall 
laminate  stiffness. 

Aboudi  (1987)  has  looked  at  a  unit  cell  representing  a  laminate  with  parallel  cracks  and 
expanded  the  displacements  using  Legendre  polynomials.  The  degraded  stiffness. was  calculated 
from  the  elastic  energy  stored  in  the  cracked  body.  The  disadvantage  of  this  method  is  that  higher 
order  Legendre  polynomials  are  required  for  good  accuracy  (Herakovich,  et  al.,  1988). 

The  principle  of  minimum  complementary  energy  was  used  to  explicitly  calculate  ply 
stresses  and  the  degraded  stiffness  by  Hashin  (1985,  1986,  1987).  This  method  is  very  accurate 
for  laminates  with  ply  stacking  sequence  of  the  form  [0P,  90q]s.  However  it  is  very  difficult  to 
extend  this  method  to  other  laminate  configurations,  because  of  the  cumbersome  nature  of  the 
complementary  strain  energy  associated  with  assumed  stress  functions  and  traction  boundary 
conditions. 

Others  have  predicted  the  degraded  stiffness  using  the  internal  state  variable  (ISV) 
concept.  This  was  originally  proposed  by  Vakulenko  and  Kachanov  (1971).  Allen  et  al. 
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(1987a,b)  have  used  this  scheme  specifically  for  transverse  matrix  cracks.  A  simple  method  for 
determining  the  ISV  was  developed  by  Lee  and  Allen  (1989)  for  predicting  the  lower  bounds  of 
degraded  axial  and  shear  moduli.  This  method  was  extended  to  include  damage  evolution  (Allen 
et  al.,  1991,  1997),  curved  matrix  cracks  (Allen  et  al.,  1988)  and  estimate  the  upper  bounds  on  the 
degraded  stiffness  (Allen  and  Lee,  1990;  Lee,  et  al.,  1991).  This  work  was  paralleled  by  Talreja 
(1985,  1986,  1989, 1992)  who  used  a  similar  approach  in  the  study  of  the  effects  and  evolution  of 

matrix  cracks. 

Another  modeling  method  using  damage  mechanics  to  describe  the  reduction  in  laminate 
stiffness  was  developed  by  Gudmundson  and  Ostlund  (1992)  and  improved  by  Gudmundson  and 
Zang  (1993).  In  this  method  the  reduction  in  stiffness  is  calculated  based  closed  form  fracture 
mechanics  solutions  for  uniform  distributions  of  cracks.  In  order  to  apply  closed  form  solutions  it 
is  assumed  that  the  cracks  are  in  an  infinite  homogeneous  medium.  This  method  gives  reasonable 
results,  but  the  accuracy  is  dependent  upon  the  specific  laminate  material  and  lay-up.  This 
method  was  extended  to  address  the  effect  on  bending  stiffness  by  Adolfsson  and  Gudmundson 
(1997).  This  is  noteworthy  since  most  models  are  applicable  only  to  pure  extension,  and  only  a 
few  address  the  changes  in  matrix  crack  behavior  under  bending  loads  (Makins  and  Adali,  1991; 
Smith  and  Ogin,  1999).  This  method  can  also  estimate  the  reduction  in  transverse  shear  stiffness. 

Less  work  has  been  reported  on  applying  the  available  models  to  actual  structures  and 
addressing  the  complications  created  by  crack  closure.  Kim,  Atluri  and  Loewy  (1998)  have 
looked  at  the  effects  of  matrix  cracking  on  flutter  for  composite  plates.  Crack  closure  is  of 
particular  concern  for  plates  under  bending  where  cracks  will  open  and  close  defending  upon  the 
deformed  shape.  Under  compression  the  matrix  cracks  will  close,  creating  a  change  in  laminate 
stiffness.  This  essentially  results  in  a  bimodular  material  with  different  stiffness  for  tension  and 
compression.  Bimodularity  is  a  phenomenon  not  exclusive  to  matrix  cracking.  Many  materials 
such  as  kevlar-rubber  composites  exhibit  this  behavior  and  the  vibration  of  bimodular  plates  has 
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been  well  studied  (Bert,  et  al„  1981;  Reddy,  1982;  Wu,  Jing  and  Cheng,  1989;  Tseng  and  Bai, 
1993).  This  work  makes  use  of  some  of  these  bimodular  methods  in  the  study  of  plates  with 
matrix  cracks. 

1.4.  Optimization 

Once  one  is  able  to  model  a  system,  such  as  an  adaptive  structural  system,  the  next 
logical  consideration  is  how  to  optimize  that  system  for  maximum  performance.  Optimization  of 
adaptive  structures  is  a  popular  area  of  study,  but  only  a  limited  amount  of  work  has  been  done  on 
simultaneously  optimizing  both  structural  and  electrical  aspects  of  a  system.  A  common 
application  of  the  electrical  interaction  with  the  structural  deformation  is  in  the  design  of  passive 
electrical  damping  systems.  The  work  by  Hagood  and  Flotow  (1991)  and  Wu  (1996,  1998)  has 
addressed  optimal  electrical  parameters  for  passive  damping  circuits.  Research  by  Tsai  and 
Wang  (1996)  has  made  use  of  optimization  techniques  to  choose  electrical  parameters  in  order  to 
maximize  control  authority. 

However,  in  all  of  these  works,  although  significant  effort  was  made  in  addressing  the 
electrical  aspects  of  the  damping  systems,  very  simple  structural  models  were  used  that  do  not 
take  into  account  the  complex  state  of  strain  that  may  exist  in  the  piezoelectric  patches.  In 
addition,  the  placement  of  the  piezoelectric  patch  has  not  been  considered  concurrently  with  the 
design  of  the  electrical  system.  Much  of  the  work  available  in  the  literature  has  addressed  the 
vibration  of  cantilevered  beams  and  plates  where  a  piezoelectric  patch  located  near  the  root  can 
be  effective  in  controlling  all  of  the  lower  order  modes.  Structures  with  other  boundary 
conditions,  where  optimal  location  varies  for  each  mode  and  are  not  intuitively  obvious  have  not 


been  examined  in  detail. 
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In  this  work  it  is  demonstrated  how  multidisciplinary  optimization  (MDO)  techniques  can 
be  utilized  to  optimize  both  structural  and  electrical  aspects  of  an  adaptive  structural  system.  To 
simultaneously  optimize  multiple  performance  requirements,  the  Kreisselmeier-Steinhauser  (K-S) 
function  approach  (Wren,  1989;  Sethi  and  Striz,  1997)  is  used.  The  K-S  technique  is  a 
multiobjective  optimization  procedure  that  combines  all  the  objective  functions  and  the 
constraints  to  form  a  single  unconstrained  composite  function  to  be  minimized.  Then  an 
unconstrained  solver  is  used  to  locate  the  minimum  of  the  composite  function.  The  advantage  of 
this  method  is  that  it  does  not  rely  on  arbitrary  weight  factors  to  combine  multiple  objective 
functions  although  they  can  be  used  in  cases  where  a  designer  wishes  to  emphasize  specific 

design  criteria  (Rajadas,  Jury  and  Chattopadhyay,  2000). 

Optimization  of  an  integrated  structural  and  electrical  system  involves  both  discrete  and 
continuous  variables.  Gradient  based  methods  are  generally  ineffective  in  the  optimization  of 
discrete  variables  and  nongradient  based  techniques  such  as  genetic  algorithms  (GA)  and 
simulated  annealing  (SA)  can  be  computationally  very  expensive  (Belegundu  and  Chandrupatla, 
1999).  A  hybrid  method  (Seeley,  Chattopadhyay  and  Brei,  1996)  has  been  developed  that  uses  a 
discrete  search  combined  with  a  gradient  based  technique  for  the  continuous  design  variables. 
The  use  of  this  technique  provides  significant  improvement  in  computational  efficiency  over 
traditional  discrete  searches  by  using  a  combinatorial  search  algorithm  that  allows  the  use  of 
gradients  for  the  continuous  variables,  while  using  a  discrete  search  technique  for  the  discrete 
design  variables.  The  optimization  method  used  in  this  work  is  very  similar  to  the  hybrid  method 
of  Seeley,  with  the  major  difference  being  that  the  discrete  search  and  the  gradient  based  search 
are  more  independent  due  to  the  nature  of  this  particular  problem.  Since,  the  present  model 
assumes  that  only  the  electrical  components  to  be  continuous  parameters  the  gradient  based 
search  essentially  amounts  to  finding  the  optimal  set  of  electrical  components  for  a  give  structural 

configuration. 


2.  Objectives 


The  overall  objective  of  this  research  is  to  develop  a  comprehensive  model  for  accurately 
and  efficiently  modeling  smart  composite  structures  for  practical  aerospace  applications  under  a 
variety  of  loads  and  environmental  conditions.  First  a  new,  efficient  method  for  modeling  smart 
structures  with  piezoelectric  devices  is  developed.  The  coupled  model  simultaneously  solves  for 
the  mechanical  and  electrical  response  of  the  system  using  mechanical  displacements  and 
electrical  displacements.  The  developed  theory  utilizes  a  refined  higher  order  displacement  field 
that  accurately  captures  the  transverse  shear  deformation  in  moderately  thick  laminates.  The 
model  is  validated  using  a  combination  of  experimental  data  and  solutions  to  known  problems. 

The  model  is  then  extended  to  include  internal  damage  in  the  form  of  delamination  and 
matrix  cracking.  When  delamination  is  present  the  sublaminates  are  modeled  as  individual  plates 
and  continuity  reinforced  at  the  interfaces.  Matrix  cracking  is  modeled  as  a  reduction  in  laminate 
stiffness  using  parameters  determined  using  finite  element  analysis  of  a  representative  crack.  The 
developed  approach  is  then  correlated  with  experimental  results.  The  effect  of  damage  on  the 
response  of  a  smart  structural  system  is  then  examined. 

Finally,  the  simultaneous  optimization  of  both  mechanical  and  electrical  parameters  in  an 
adaptive  structural  system  is  studied.  The  objective  of  this  work  is  to  demonstrate  how 
multidisciplinary  optimization  (MDO)  techniques  can  be  utilized  to  optimize  both  structural  and 
electrical  aspects  of  an  adaptive  structural  system.  To  simultaneously  optimize  multiple 
performance  requirements,  the  Kreisselmeier-Steinhauser  (K-S)  function  approach  is  used. 
Optimization  of  piezoelectric  actuator  placement  and  electrical  circuitry  is  performed  on  passive 


electrical  damping  systems. 
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The  following  are  the  principle  objectives  for  this  research. 

1 .  Development  of  an  efficient  coupled  model  for  simultaneous  modeling  of  mechanical  and 
electrical  response  in  an  adaptive  structural  system. 

2.  Development  of  a  unified  method  for  modeling  the  effects  of  delamination  on  the  static  and 
dynamic  response  of  smart  composite  structures. 

3.  Development  of  a  model  for  predicting  the  reduction  in  plate  stiffness  caused  by  transverse 
matrix  cracking,  and  a  method  for  describing  its  effects  on  the  structural  response  of  smart 
composites. 

4.  Analysis  of  the  effects  of  combined  damage  on  smart  structural  response  (both  static  and 
dynamic)  using  the  developed  procedure. 

5.  Demonstrate  the  use  of  optimization  techniques  in  the  simultaneous  optimization  of 
mechanical  and  electrical  components  in  a  smart  structural  system. 


3.  Piezoelectric  adaptive  structures  modeling 


In  this  chapter  a  coupled  piezoelectric-mechanical  model  is  developed  in  order  to 
accurately  and  efficiently  model  both  the  mechanical  and  electric  response  of  an  adaptive 
structural  system.  Unlike  other  techniques  that  solve  only  one  of  the  coupled  piezoelectric 
equations,  the  developed  model  simultaneously  solves  both  of  the  coupled  equations.  This 
ensures  that  the  model  accurately  describes  the  transfer  of  energy  between  mechanical  and 
electrical  and  that  the  model  will  be  applicable  regardless  of  whether  the  piezoelectric  device  is 
being  used  as  a  sensor  or  actuator.  Unlike  other  models,  the  developed  theory  ensures 
conservation  of  charge  and  accurately  predicts  the  electrical  response  of  the  system.  Since  the 
electrical  response  of  the  system  is  considered,  any  system  using  piezoelectric  materials  can  be 
modeled  other  than  just  the  tradition  sensor  and  actuator  applications.  These  systems  include 
passive  electrical  damping  circuits  and  self-sensing  actuators. 

A  less  frequently  used  form  of  the  coupled  constitutive  relations  is  shown  to  be 
advantageous  for  analysis  for  adaptive  structures  and  is  used  to  develop  a  linear  model.  A  refined 
higher  order  laminate  theory  is  used  to  model  the  mechanical  displacement  field,  thus  making  the 
model  capable  of  capturing  the  effects  of  transverse  shear  in  moderately  thick  laminate 
constructions.  Additional  equations  are  then  included  to  model  the  attached  electrical  circuitry. 
The  resulting  model  is  verified  using  experimental  data,  and  the  importance  of  solving  the 
coupled  system  is  illustrated  in  the  transient  analysis  of  adaptive  structural  systems. 


Fig.  3.1.  Plate  geometry  and  coordinate  system. 


If 
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3.1.  Linear  coupled  piezoelectric-mechanical  theory 

The  model  developed  in  this  section  is  based  on  solving  the  response  for  a  laminated 
composite  plate  with  embedded  or  surface  bonded  actuators.  A  two-dimension  plate  theory  is 
used  for  computational  efficiency.  The  plate  geometry  is  illustrated  in  Figure  3.1.  The  global  x- 
and  y-  coordinates  are  assumed  to  be  oriented  in  the  plane  of  the  plate,  and  the  z-  axis  is 

perpendicular  to  the  plane  of  the  plate.  The  plate  is  a  composite  laminate  composed  of 

orthotropic  laminae.  The  material  1-  and  2-  axes  are  assumed  to  be  parallel  to  the  plane  of  the 
plate,  and  the  material  3-  axis  is  oriented  parallel  with  the  global  z-  axis. 

3.1.1.  Constitutive  relations 

The  construction  of  a  model  for  smart  composite  laminates  starts  with  the  formulation  of 
the  constitutive  relations  (IEEE,  1987).  Traditionally  these  are  expressed  as  a  function  of  the 
components  of  strain  (S9)  and  electric  field  (£,)  as  follows 

Ty  =  cfjklSkl  -  ekijEk  (3-D 

'  Di=eildSld  +  zfkEk  (3-2) 

where  T„  is  the  stress  tensor,  D,  is  the  electric  displacement,  and  cEijkh  ekiJ  and  %sik  are  elastic, 

piezoelectric  and  dielectric  constants  for  the  material.  The  superscript,  E,  in  the  elastic  constant 
reflects  that  the  elastic  constants  used  are  those  measured  in  the  absence  of  any  electric  field 
within  the  piezoelectric  material.  Likewise,  the  superscript,  5,  in  the  dielectric  constant 
represents  those  constants  measured  while  the  piezoelectric  material  is  subjected  to  zero  strain. 
The  electric  field  is  generally  thought  of  as  being  derived  from  the  gradient  of  a  scalar  electric 


potential,  <f>  as  follows 
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t 


E,=-+«  (3-3) 

When  piezoelectric  materials  are  used  in  adaptive  structural  systems  they  are  usually  in 
the  form  of  relatively  thin  layers  that  are  either  embedded  in  or  bonded  to  the  surface  of  the  host 
structure.  The  upper  and  lower  surfaces  are  usually  each  covered  with  an  electrically  conducting 
electrode  in  order  to  apply  or  sense  changes  in  voltage,  while  the  sides  generally  remain 
unelectroded.  Based  on  this  typical  configuration  structural  modeling  has  traditionally  uncoupled 
the  constitutive  relations  by  assuming  that  the  electric  field  through  the  thickness  of  the 
piezoelectric  device  is  directly  proportional  to  the  voltage  on  the  electrodes  divided  by  the 
material  thickness.  Having  made  this  assumption,  Eq.  (3.1)  can  directly  solved  for  in  order  to 
predict  the  response  of  the  structure  in  the  presence  of  an  applied  voltage.  This  has  often  been 
referred  to  as  the  converse  effect.  In  a  similar  manner,  Eq.  (3.2),  known  as  the  direct  effect,  can 
be  directly  solved  for  in  order  to  give  the  charge  flow  into  or  out  of  the  piezoelectric  device  in  the 
presence  of  an  applied  strain.  This  allows  the  piezoelectric  sensor  to  give  an  estimate  of  the  local 
strain  in  the  host  structure. 

While  the  uncoupled  approximation  gives  convenient  estimates  concerning  the  response 
of  an  adaptive  structural  system,  it  ignores  some  key  physical  principles  of  piezoelectric 
materials.  First,  the  uncoupled  approach  ignores  the  nature  in  which  energy  is  transformed  from 
electric  to  mechanical  and  vice  versa.  When  a  voltage  is  applied  to  the  electrodes  of  an  actuator, 
an  electric  field  is  induced  which  in  turns  creates  strain  in  the  host  structure.  However,  the  strain 
in  the  piezoelectric  material  also  influences  the  electric  field,  as  shown  in  Eq.  (3.2).  Thus,  the 
electric  field  in  the  actuator  is  dependent  on  the  stiffness  of  the  host  structure,  which  is  not 
captured  by  the  uncoupled  assumptions. 

Another  aspect  not  addressed  by  the  uncoupled  approach  is  conservation  of  charge. 
Physics  requires  that  the  electric  displacement  vector,  D„  satisfy  the  electric  equation  for  an 


insulator 
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This  equation  is  simply  a  mathematical  representation  of  the  conservation  of  charge  within  the 
material.  For  the  case  of  the  piezoelectric  device  with  electrodes  on  the  upper  and  lower  surfaces, 
it  implies  that  the  electric  displacement  be  constant  through  the  thickness  of  the  device.  This  is 
very  important,  in  that  Eq.  (3.2)  now  shows  that  when  the  strain  varies  through  the  thickness, 
such  as  during  plate  bending,  then  the  electric  field  must  also  vary  through  the  thickness  in  order 
to  maintain  constant  electric  displacement.  This  also  implies  that  the  uncoupled  modeling 
approach  results  in  differing  amounts  of  charge  on  the  upper  and  lower  electrodes,  which  is  a 
violation  of  the  conservation  of  charge  principle.  This  has  been  resolved  by  some  authors  by 
solving  the  coupled  equations,  but  making  the  electric  potential,  and  in  turn  the  electric  field, 
higher  order  functions  of  the  through  the  thickness  coordinate  to  match  the  displacement  and 
strain  fields  in  the  structure.  However,  such  an  approach  can  lead  to  numerous  additional  degrees 
of  freedom  to  describe  the  electric  potential.  Another  drawback  of  such  an  approach  is  that  the 
resulting  system  matrices  in  finite  element  implementation  are  not  symmetric.  This  results  in  a 
sizable  increase  in  the  computational  effort  required  to  solve  the  system  of  equations  when 
compared  to  the  original  host  structure  whose  finite  element  equations  result  in  a  symmetric 
positive  definite  system. 

Equations  (3.1)  and  (3.2)  are  not  the  only  way  to  describe  a  piezoelectric  material.  The 
following  forms  are  all  equally  valid 


Ty  -  cDijkiSkl  hkijDk 

(3.5a) 

Ei=-hmSki+ps*Dk 

(3.5b) 

Sy  =  sE  ijkiTfj  +  dkij  Ek 

(3.6a) 

Di  =  diklTkl  +  x  ik E k 

(3.6b) 

t 
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=  sDijidTkl  +gUjDk 

(3.7a) 

-SuTu  +PT*Dk 

(3.7b) 

While  each  form  of  the  coupled  piezoelectric  equations  has  advantages  in  certain  cases,  the  form 
using  the  mechanical  strain  and  the  electric  displacement,  Eqs.  (3.5a,  3.5b),  are  best  suited  to 
adaptive  structural  modeling.  The  reason  for  this  is  that  standard  finite  element  solutions  for  the 
mechanical  response  of  structures  utilize  displacements  to  provide  strain,  and  specifying  electric 
displacement  as  an  in-plane  variable  inherently  satisfies  conservation  of  charge.  Also,  the 
resulting  system  of  equations  is  symmetric  allowing  for  improved  computational  efficiency.  It  is 
also  worth  noting  that  when  the  piezoelectric  device  is  used  as  an  actuator,  the  applied  voltage  is 
introduced  in  a  manner  very  similar  to  how  external  forces  and  pressures  are  applied  to  the 
structure. 

The  constitutive  relations  can  be  expressed  in  matrix  form  as 


cr  =  cDf-hTD 


(3.8a) 


E  =  -hf  +  psD 


(3.8b) 


where  the  stress  and  strain  have  been  expressed  as  follows 
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a  =  < 


(3.9) 
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a  =  4 
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The  nature  of  the  matrices  formed  by  the  elastic,  piezoelectric  and  dielectric  constants  is 
very  much  dependent  on  nature  and  orientation  of  the  piezoelectric  crystal.  Structural  analysis  is 
simplified  by  the  assuming  that  the  piezoelectric  device  possesses  a  minimum  amount  of 
symmetry,  making  it  at  least  orthorhombic.  This  means  that  the  material  is  assumed  to  have  three 
mutually  perpendicular  crystal  axes  with  the  polar  axis  oriented  in  the  out-of-plane  direction. 
This  assumption  is  not  very  restrictive,  in  that  almost  all  piezoelectric  patches  used  in  adaptive 
structural  systems  fall  in  this  category.  It  also  assumed  that  the  entire  upper  and  lower  surfaces 
are  covered  by  perfectly  conducting  electrodes  and  that  the  sides  are  unelectroded.  This  is  the 
typical  configuration  for  piezoelectric  devices  and  results  in  the  requirement  that  the  in-plane 
components  of  the  electric  displacement  be  uniformly  zero  at  all  locations. 

D,  =  D2  =  0  (3.11) 
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(3.10) 


3. 1.2.  Use  of  material  properties 

The  components  of  the  elastic,  piezoelectric  and  dielectric  matrices  must  be  derived  from 
experimental  data.  Although  the  matrices  in  Eqs.  (3.8a,b)  can  be  measured  directly,  the  values 
usually  measured  and  provided  by  manufacturers  are  the  independent  constants  that  make  up  the 
components  the  form  of  the  constitutive  relations  from  Eqs.  (3.6a,b).  In  matrix  form  these 
equations  are  expressed  as 


e  =  sEcr  +  dTE 


(3.12a) 
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D  =  do-  +  xTE  (3.12b) 

As  a  result  of  the  orthorhombic  assumption  the  elastic,  piezoelectric  and  dielectric  matrices  take 
the  forms 


(3.13) 


0  0  0  0  dt5  0 

d  =  0  0  0  d2A  0  0  (3.14) 

_J31  dn  d3l  0  0  0_ 

>11  0  0  ’ 

X=  0  Z22  0  (3.15) 

.  0  0  ^33. 

and  all  of  the  unknown  coefficients  are  those  typically  determined  from  experimental  data.  Thus, 
in  order  to  determine  the  matrix  coefficients  in  the  form  of  the  equations  expressed  in  Eqs. 
(3.8a,b)  a  transformation  is  required.  This  aspect  of  piezoelectric  materials  is  often  not  discussed 
by  authors,  but  is  extremely  important  in  order  develop  and  accurate  set  of  system  matrices.  The 
transformation  is  complicated  by  the  fact  that  most  two-dimensional  structural  models  make  the 
assumption  of  plane  stress.  The  refined  higher  order  laminate  theory  discussed  later  in  this  work 


also  results  in  the  out-of-plane  stress  being  zero,  or 

cr3  =  0 


(3.16) 
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It  is  very  important  to  properly  address  this  since  the  piezoelectric  coefficient,  d33,  is  typically 
very  large.  Ignoring  the  effect  of  Eq.  (3.16)  on  the  piezoelectric  behavior  will  result  in  material 
coefficients  that  are  very  inaccurate.  Correct  transformation  of  results  in  the  following  equations 
for  determining  the  coefficients  for  the  matrices  in  Eqs.  (3.8a, b)  from  the  typical  matenal 

properties  available. 


„  Yfixi-YfAu) 

C"=  A 

(3.17a) 

D  Y‘(xl,-Y‘dl ,) 

22  A 

(3.17b) 

n  Y,E(Yfd,xd  »+vfa») 

Cn= - X 

(3.17c) 

n  Y^Yfd^dn+v^xl)  _  » 
c21- - A 

(3. 17d) 

D  ^23/^22 

*"44  ~  T  j2  r,E 

Xl2  ~d24G23 

(3.17e) 

D  G*Xu 

Cs5  XTn-d\filx 

(3.17f) 

c6°6  =  G* 

(3  - 1 7g) 

(3.18a) 

P22  “  yT  _ d 2  GE 

Xl2  °’24KJ23 

(3.18b) 

pS  1  ^12^21 

A 

(3.18c) 

r‘k%,+dj 

n,,  — 

(3.19a) 

h32  = 


23 

^2  (^11^31  +  ^32  ) 

(3.19b) 

A 

^15^31 

(3.19c) 

■  Xu  -4^1 

^24^23 
.J  j2  nE 

(3.19d) 

X27  aiA{Jn 


where 


a  =  xl,{ i  -  « )- r‘d2» -r/4  -  vfaAtf  -  y&M"  <3-20> 

All  other  coefficients  are  zero  or  are  not  used  by  the  refined  higher  order  laminate  theory. 


3.1.3.  Equations  of  motion 

The  conservation  of  energy  within  piezoelectric  material  is  described  by  the  electric 
enthalpy,  which  is  typically  expressed  in  terms  of  strain  and  electric  field  as 


Note  that  this  is  different  from  the  total  potential  energy,  which  is 


1  £  l-T  (-T  .  1  .S 


(3.21) 


U=-cLStSu+-z;£iEj 


(3.22) 


~  2  ,JKl  v  ki  2  J  J 

Since  the  formulation  used  in  this  work  is  based  on  strain  and  electric  displacement,  an  alternate 
form  of  the  electric  enthalpy  must  be  used 

(3.23) 


or  in  matrix  form 


R(s,  D)  =  -  f  TcDf  -  Drhs  +  ^dtpsd 
2  2 


(3.24) 


* 
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The  equations  of  motion  can  be  formulated  using  a  variational  approach  and  Hamilton’s 
Principle  in  a  maimer  similar  to  that  proposed  by  Tiersten  (1967).  The  variational  principle 
between  times  to  and  t,  for  the  piezoelectric  body  of  volume  Fean  be  written  as  follows 


sn=o=nv 


fi  . T • 

[2**  u 


•<5H(e  ,D) 


dVdt  +  £  SWdt 


(3.25) 


where  the  first  term  represents  the  kinetic  energy,  the  second  term  the  electric  enthalpy  and  5W  is 
the  total  virtual  work  done  on  the  structure.  The  terms  u  and  p  correspond  to  the  mechanical 
displacement  and  density,  respectively.  The  work  done  by  body  forces  (fB),  surface  tractions  (fs) 
and  electrical  potential  ($  applied  to  the  surface  of  the  piezoelectric  material  can  be  expressed  by 

<5W  =  j^SuTfBdV  +  £  5uTfsdS  +  £<SDT<zW.S'  (3-26) 


Equations  (3.24-26)  provide  the  equations  of  motion  for  the  piezoelectric  body.  To  solve  them, 
the  assumptions  about  the  properties  and  orientation  of  the  piezoelectric  device  described  in  the 
previous  sections  are  utilized  to  simplify  the  system. 


3. 1.4.  Refined  higher  order  laminate  theory 

A  refined  higher  order  laminate  theory,  originally  developed  by  Reddy  (1984,  1997),  is 
used  to  model  the  mechanical  displacement  field.  The  refined  higher  order  theory  assumes  a 
parabolic  distribution  of  transverse  shear  strain,  thus  providing  accurate  estimation  of  transverse 
shear  stresses  for  moderately  thick  constructions  with  little  increase  in  computational  effort.  This 
makes  the  higher  order  theory  more  accurate  than  the  classical  plate  theory,  which  neglects  the 
presence  of  transverse  shear,  and  more  convenient  than  the  first  order  shear  deformation  theory, 

which  requires  a  shear  correction  factor  to  be  know. 

The  laminate  is  assumed  to  be  a  plate  structure  composed  of  an  arbitrary  number  of 
orthotropic  lamina  arranged  with  varying  orientations.  The  coordinate  system  for  the  plate  is 
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taken  to  be  with  the  x-y  plane  parallel  to  the  plane  of  the  plate  and  the  z  coordinate  normal  to  the 
plane  of  the  plate  measured  from  the  center  of  thickness.  The  theory  starts  with  a  general,  third 


order  displacement  field  as  follows 


u,  =  K  +  z^,+zVJ  2  +2V*3 


U2  =V  +  Z^,+ZV,2  +ZVj,3 


u3  =  w 


(3.27a) 


(3.27b) 


(3.27c) 


The  variable  z  represents  the  location  with  respect  to  the  midplane  of  the  plate,  and  h  is  the  total 
plate  thickness.  This  displacement  field  is  then  simplified  by  imposing  the  stress  free  boundary 


conditions  on  the  free  surfaces. 


Ml  Mi  n 


(3.28) 


Since  the  laminate  is  orthotropic,  this  implies  that  the  transverse  shear  strains  are  zero  on  the 


upper  and  lower  surfaces  of  the  laminate. 


\  2)  dz{  2)  dy  {  2) 
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(3.29a) 


(3.29b) 


Thus,  substituting  Eqs.  (3.27a-c)  into  Eqs.  (3.29a,b)  yields  the  following 

=  Vyi  =  0 


(3.30a) 


-4  (  dw 


(3.30b) 


-4  f  dw 

y y*  ~  31,2  ^ y 1  +  pht 


(3.30c) 


In  order  to  simplify  the  resulting  equation  the  following  change  of  variables  is  made 
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dw 

(3.31a) 

<i|<& 

+ 

II 

(3.31b) 

The  refined  displacement  field  now  takes  the  following  form 


cHv'i  4  zi 


U'a“  +  V'"&J"3Fr‘ 

f  4  z3 

“2=v  +  z  n-- 


(3.32a) 

(3.32b) 

(3.32c) 


where  u  ,  v,  and  w  are  the  displacements  of  the  midplane  and  y/x  and  i//y  are  the  rotations  of  the 
normal  at  z=0  about  the  -y  and  -x  axes,  respectively.  Note  that  u,  v,  w,  y/x  and  y/y  are  all  functions 
of  the  x  and  y  coordinates  only.  The  coupled  piezoelectric-mechanical  theory  can  easily  be 
developed  using  other  plate  theories,  but  the  chosen  plate  model  affects  the  constraints  imposed 
on  the  delamination  boundaries  developed  later  in  this  paper.  The  y/x  and  \f/y  values  can  be 
thought  of  as  quantifying  the  amount  of  transverse  shear  present  in  the  laminate,  and  if  they  are 

assumed  to  be  zero  the  classical  plate  theory  is  obtained. 

The  linear  form  of  the  strains  are  then  determined  by  differentiating  the  displacements, 

resulting  in  the  following  relations 


du]  _du  dy d2w)  4 z3  dy^ 
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(3.33d) 

(3.33e) 

(3.33f) 


3. 1.5.  Finite  element  formulation 

To  solve  the  equations  of  motion  for  the  coupled  system  a  finite  element  solution  method 
is  used.  This  method  is  chosen  in  order  to  allow  relatively  easy  modeling  of  a  variety  of 
structural  geometries.  Since  the  coupled  approach  that  is  being  used  solves  for  both  mechanical 
and  electric  aspects  of  the  coupled  system,  the  finite  element  method  used  to  solve  these 
equations  must  address  both  the  mechanical  and  electric  variables  in  order  to  be  effective. 


3. 1.5.1  Mechanical  relations 


By  using  the  refined  higher  order  displacement  field,  the  equations  of  motion  can  be 
simplified  by  the  reduction  of  the  mechanical  displacements  and  strains  to  the  variables  u,  v,  w, 
y/x  and  \f/y.  The  displacement  vector,  u,  can  be  written  as  follows 


u  =  < 


u, 


u, 


*3J 


r  =  ®uUul 


(3.34) 


where 


«ul  = 


U  ¥x  V  ¥y 
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dw  dw 
dx  dy 


(3.35) 
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(3.36) 


Since  the  laminae  in  the  plate  are  all  assumed  to  be  orthotropic  the  m-plane  strains  and  the 
transverse  shear  strains  can  be  decoupled  and  treated  independently  from  on  another.  As  a  result 
of  the  plane  stress  assumption,  the  out-of-plane  strain,  s3,  is  ignored.  This  simplifies  the 
assembly  process  slightly  and  gives  two  separate  strain  vectors.  The  in-plane  linear  strain  vector, 
g|L,  and  the  out-of-plane  strain  vector,  Eql>  are  defined  as  follows 


e„  =  =B„b 


(3.37) 


£ol  ~ 


B01uul 


(3.38) 


where 


&  &  ^  aSv  (3  39) 

~dx  dy  dx  dy  dx  8y  dx  dy  dx 2  dy2  dxdy  \ 
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As  stated  earlier  the  in-plane  components  of  the  electric  displacement  are  known  to  be 
zero  at  all  locations,  thus  only  the  out-of-plane  electric  displacement,  D3,  is  left  as  an  unknown. 
Substitution  of  Eqs.  (3.34,37,38)  into  the  equations  of  motion,  Eqs.  (3.24-26),  yields  the 
following  system  of  equations 

-  .puXpBA,  +<5uLBiilc“B,lU.!  +  ouI1BJjlc"B0lu»,L/a  = 
J.lL-<Si:XLhTD,-a>,hB,Lu.J+a),p“3D,  J  (342) 

( li  *1,B  IhdV  +  ( {<S>Xfs  +  +  &lPTJr  V 

Since  all  unknowns  are  functions  of  the  in-plane  coordinates,  x  and  y,  only,  integration  through 
the  thickness  is  performed.  It  should  be  noted  that  the  integration  between  the  lower  and  lower 
surfaces,  zL  and  zv,  actually  must  involve  a  summation  of  integrations  over  each  individual 
lamina,  since  the  material  properties  can  not  be  assumed  to  be  constant  through  the  thickness  of 
the  plate.  First,  the  inertia  term  can  be  reduced  to 

£1'b>B,*  =  L>„  (3.43) 


where 


o  II, 
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(3.44) 
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(3-45) 
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1  0  0  0  0  0  0' 

0  0  1  0  0  0  0 

0  0  0  0  1  0  0 

0  1  0  0  0  -1  0 

0  0  0  1  0  0  -1  (3.46) 
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Next,  the  body  force  can  be  reduced  by  assuming  that  the  local  force  is  proportional  the  local 

density,  as  is  generally  the  case  for  inertial  loading. 

f.  =  A  <3-47) 

Thus,  the  body  force  term  becomes 

fLBXdz  =  LlUB  (3-48) 

where 

X 

I  =  II,  (3-49) 

X. 

Through-the-thickness  integration  of  the  strain  potential  energy  terms  yield 

£UB7lc|,Bil&  =  XCjLil  (3.50) 

Bolco®ol^z  =  LolCoLol 

where 
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ci  ct  ct_ 


(3.51) 
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Through-the-thickness  integration  of  the  piezoelectric  coupling  terms  yield 

£ub>td3Jz  =  lTlptd3 


where 


p=[p.  p,  pj 

P(  =  £  h z'dz 
h  =  [/t31  hn  0] 
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(3.52) 

(3.53) 

(3.54) 
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Through-the-thickness  integration  of  the  electric  potential  energy  terms  yields 

£"D,pf,DJ*  =  D,D  (361) 

where  t„,  is  the  thickness  of  the  piezoelectric  layer.  The  displacement  vectors  in  Eqs.  (3.35)  and 
(3.39)  can  be  expressed  using  a  finite  element  interpolation  from  the  nodal  vector,  ue. 

Uu,=Nul(^^)ue  (3-62) 

u„2  =Nu2(*,)/)Ue  (3.63) 

The  finite  element  scheme  developed  in  this  work  uses  linear  interpolation  of  the  variables  u,  v, 
yrx  and  y/y,  and  a  Hermit  cubic  polynomial  function  for  the  out  of  plane  displacement,  w.  This 
results  in  seven  mechanical  degrees-of-freedom  per  node,  u,  v,  if/x,  y/y,  w,  dw/dx,  and  dw/dy. 


3. 1.5. 2  Electric  displacement 

Now  that  the  mechanical  displacements  have  been  reduced  via  finite  elements  to  nodal 
degrees  of  freedom,  it  is  now  necessary  to  address  the  distribution  of  electric  displacement.  Since 
electrodes  cover  the  entire  in-plane  surfaces  of  the  piezoelectric  layer,  the  electric  displacement 
can  be  described  as  the  charge  per  unit  area  on  the  electrodes. 

n  =^i  (3-64) 

3  dA 

where  q  is  the  amount  of  positive  charge  on  the  lower  electrode  and  A  is  the  surface  area  of  the 
PZT  electrode.  There  are  two  possible  methods  for  treating  the  distribution  of  electric 
displacement  across  the  piezoelectric  layer.  The  most  general  method  is  to  use  a  finite  element 
interpolation,  and  relate  the  distribution  to  the  electric  displacement  at  a  given  set  of  nodes.  The 
electric  displacement  then  becomes 


D3  =Nq(x,y)De 


(3.65) 
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where  De,  is  the  nodal  electrical  displacement  and  Nq(x,>0  is  the  interpolation  matrix.  Integrating 
over  the  area  of  the  plate  using  Gauss  quadrature  on  each  of  the  finite  elements  yields  gives  the 
governing  equations  for  the  coupled  theory. 

<5n„  =<5uJ  f(-M„U,-KraU«-K„DD,+Fu)rf(  =  0  (3.66a) 

<S1D  =  f  (-K„uii,  -KDDDe  +F B)dt  =  0  (3.66b) 

Jto 

where  Mu  is  the  structural  mass  matrix,  the  matrix  Kuu  is  the  mechanical  stiffness  matrix,  Kdd  is 
the  electrical  stiffness  matrix,  and  Kud  and  KDu  are  the  stiffness  matrices  due  to  piezoelectric- 
mechanical  coupling.  The  vectors  Fu  and  FD  are  the  force  vectors  due  to  mechanical  and 
electrical  loading.  These  matrices  and  vectors  are  defined  as  follows 


ngp  ngp  _ 

Mu  =  J|NIiLIILuNui 

/=1  j=\ 

(3.67) 

Kuu  =  yyw,.w/|jl{NLL^C1L1LNU2  +Nu1LolC0LolNu,} 

<= i  M 

(3.68a) 

ngp  ngp 
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ngp  ngp 
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where  ngp  is  the  number  of  gauss  points  used  and  w,-  and  wj  are  the  gauss  point  weights.  The 
value  Ae  represents  the  area  of  each  element  and  is  necessary  for  converting  electric  displacement 

to  total  charge. 

The  result  system  of  equations  can  be  expressed  in  matrix  form  as 


Mu  0  fiiel  Cu  0  fue]  Kuu  -^ud 

0  0  [DJ  0  0  [I)]  _^Du  ^DD 


(3.70) 


To  incorporate  structural  damping  into  the  equations  a  structural  damping  matrix  C„  is  added. 
The  nature  of  the  damping  matrix  can  be  chosen  to  meet  the  needs  of  the  user.  In  this  work 
classical  damping  has  been  used. 

A  second  method  is  possible  for  cases  when  the  strain  across  the  piezoelectric  device  can 
be  assumed  to  be  constant.  This  assumption  would  be  valid  for  small  patches  attached  to 
relatively  large  structures  and  not  near  any  stress  concentrations.  Since  the  strain  across  the 
piezoelectric  layer  is  assumed  to  be  constant,  then  the  electric  displacement  can  thereby  by 
assumed  to  be  constant.  For  this  case  Eq.  (3.64)  is  reduced  to 

D,  =  —  «—  forsmallPZTs  (3-71) 

3  a 4  a 

where  q  is  the  total  charge  on  the  electrodes.  Thus,  the  governing  equations  can  now  be 
expressed  in  terms  of  the  charge  flow  in  each  piezoelectric  device.  The  advantage  of  this  is  that 
each  piezoelectric  device  added  to  the  structure  only  adds  a  single  electrical  degree  of  freedom  to 
the  system.  Thus,  modeling  of  a  complex  smart  structural  system  can  be  accomplished  with  little 
more  computational  effort  that  that  required  for  the  original  structure.  Integrating  over  the  area  of 
the  plate  yields  gives  the  governing  equations 

<5n„  =  <5uJ  f  (-M.fi,  -Ku„u,  -K„q  +  Fu)tf  =  0  (3.72a) 


<*1,  =  V  ('  (-K„u,  -K„q  +  F,)*  =  0 


(3.72b) 
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where  q  is  the  vector  of  the  PZT  charges  and  Mu  is  the  structural  mass  matrix,  and  Cu  is  the 
structural  damping  matrix.  The  matrix  Kuu  is  the  mechanical  stiffness  matrix,  Kqq  is  the  electrical 
stiffness  matrix,  and  Kuq  and  K„u  are  the  stiffness  matrices  due  to  piezoelectric-mechanical 
coupling.  The  vectors  Fu  and  Fq  are  the  force  vectors  due  to  mechanical  and  electrical  loading. 
These  terms  are  defined  in  Eqs.  (3.67),  (3.68a),  (3.69a)  and  as  follows 


1  ngp  ngp 

k„  =-^2MjInW..p 

A  /=i  M 

(3.73a) 

1  ngp  ngp 

K,u=K;=-ixi^lJlPL-N- 

A  M  >1 

(3.73b) 

1  ngp  ngp 

k  = — yyw 
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A  i-i  >i 

(3.73c) 
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(3.74) 

Similarly,  the  governing  equations  can  be  written  in  matrix  form  as  follows 


Mu  0 
0  0 


u„ 


•  + 


Cu  0 
0  0 


u. 


K  K 

“uu  xvuq 

K  K 

^qu  X^qq 


(3.75) 


The  choice  of  which  method  to  use  depends  on  the  structural  system  being  analyzed.  If  a 
strain  gradient  across  the  piezoelectric  layer  is  present  then  large  errors  can  result  from  using  the 
second  method.  For  the  majority  of  the  results  shown  in  this  work  the  finite  element 
discretization  of  the  electric  displacement  has  been  used  to  ensure  accuracy. 


3.1.6.  Electrical  equations 

The  absence  of  any  electrical  inertia  or  damping  terms  in  Eqs.  (3.70)  and  (3.75)  is  a  result 
of  only  considering  the  mechanical  aspects  of  the  smart  structure.  When  considering  a  smart 
structure  as  a  whole,  additional  terms  must  be  added  for  electrical  components  in  the  system.  For 
a  simple  LRC  circuit,  the  variational  energy  must  include  the  following 
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<5nq  tW 

V  2  y 


+  F<% 


(3.76) 


where  L,  R  and  C  are  the  inductance,  resistance  and  capacitance,  and  V  is  the  applied  voltage. 
The  importance  of  formulating  Hamilton’s  Principle  in  terms  of  the  charge,  rather  than  electric 
field  or  potential,  now  becomes  apparent.  If  the  equations  of  motion  for  any  electrical  system 
attached  to  the  smart  structure  can  be  formulated  in  the  following  form 

Mqqe  +  Cqqe  +  Kqqe  =  Fq  (3.77) 

then  these  equations  can  be  directly  combined  with  Eq.  (3.75).  This  combination  results  in  a 
completely  coupled  electrical-mechanical  system  of  the  following  form 


where  qe  includes  not  only  the  charge  associated  with  the  PZTs,  but  also  the  electrical  system. 
Matrix  Kqq*  also  includes  the  electrical  stiffness  components  of  both  the  PZTs  and  the  electrical 
system.  If  the  electric  displacement  was  discretized  then  it  is  necessary  to  use  the  following 
relation  between  the  electric  displacement  and  the  total  PZT  charge 


?,=(£Vs)d,.  or  q  =  A„D,  (3.79) 

Combining  Eqs.  (3.70),  (3.77)  and  (3.79),  the  resulting  coupled  electrical-mechanical 
system  equations  are  obtained 

Mu  0  IJiil  rc„  0  IJu.l  IX.  K«D  "IHLJM  (3.80) 
o  aXaJId.JI0  a;c,a,J\dJ  |k„.  kdd+aJk,a,J\dJ  [fJ 

Equations  (3.78)  and  (3.80)  provide  the  means  to  solve  a  variety  of  coupled  piezoelectric- 
mechanical  problems.  The  eigenvalues  and  eigenvectors  can  be  determined  to  compute  the 
natural  frequencies  and  mode  shapes  of  the  system.  The  complex  linear  problem  can  be  solved  to 
determine  the  forced  response  of  the  system  to  a  harmonic  excitation.  Another  advantage  of  this 
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model  is  that  the  integration  of  the  structural  finite  element  model  and  electrical  components 
allows  for  design  of  passive  damping  circuits  for  complicated  structures  using  multi-variable 
optimization  methods.  This  technique  would  allow  determination  of  optimum  PZT  location, 
orientation  and  the  electrical  component  values  for  single  or  multi-mode  damping  of  the  structure. 


3. 1. 7.  Boundary  Conditions 

The  developed  system  equations  shown  in  Eqs.  (3.78)  and  (3.80)  apply  to  an 
unconstrained  plate  and  do  not  take  into  account  any  mechanical  or  electrical  constraints. 
Mechanical  constraints  are  applied  to  the  finite  element  system  in  the  conventional  manner  using 
the  elimination  approach.  Boundary  conditions  for  the  electrical  system  must  also  be  applied 
depending  on  the  nature  of  the  electrical  system  that  is  assumed. 

If  the  electrical  system  is  fully  modeled  using  Eq.  (3.78)  or  (3.80)  then  the  electrical 
system  is  fully  described  and  additional  boundary  conditions  are  generally  not  required. 
However,  it  is  also  permissible  to  not  include  the  electrical  circuitry  equations  and  instead  use  Eq. 
(3.70)  or  (3.75)  to  describe  the  system.  In  this  manner  a  direct  voltage  can  be  assumed  to  be 
applied  to  each  piezoelectric  device  or  a  fixed  amount  of  charge  flow  can  be  applied  as  a 
boundary  condition.  This  method  neglects  the  effects  of  resistance  and  inductance  in  the 
electrical  system,  but  it  is  in  most  cases  convenient  and  accurate.  The  application  of  a  fixed 
voltage  corresponds  to  the  traditional  method  of  in  which  an  actuator  is  controlled.  If  the  applied 
voltage  is  zero  then  the  PZT  is  assumed  to  be  shorted  out  and  the  resulting  charge  flow  from  the 
device  acts  as  a  sensor  signal.  The  voltage  is  applied  as  a  load  and  appears  in  the  forcing  vector 
on  the  right-hand  side.  If  the  piezoelectric  device  is  assumed  to  be  open  circuited  then  the  net 
charge  flow  from  the  device  will  be  zero.  If  the  system  is  formulated  in  terms  of  total  charge 
using  Eq.  (3.75),  then  the  fixed  value  of  charge  is  easily  applied  as  a  boundary  condition  using  the 
elimination  method.  However,  if  the  electric  displacement  has  been  discretized  using  nodal 
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electrical  displacements  with  Eq.  (3.70)  then  one  must  be  careful  to  specify  that  the  net  charge 
flow  be  zero,  not  the  individual  nodal  electrical  displacements.  The  net  charge  flow  in  Eq.  (3.79) 
is  specified  to  be  zero  creating  a  single  constraint  equation  for  each  piezoelectric  device.  This 

equation  can  be  enforced  using  the  penalty  method  creating  a  submatrix,  cpen  Aq  Aq ,  to  be  added 

to  the  system  stiffness  matrix.  Results  show  that  the  solution  is  rather  insensitive  to  the  penalty 
value,  Cpen,  and  that  values  around  102  times  the  maximum  diagonal  component  ofKuu  work  well. 

Once  appropriate  boundary  conditions  are  applied  the  resulting  system  of  equations  can 
be  solved  based  on  the  type  of  analysis  required.  For  ease  of  notation,  this  system  of  equations  is 
simplified  to  the  following  form 

ma  +  cv  +  ku  =  f  (3-81) 

It  is  important  to  note  that  the  system  vectors  u,  v,  and  a  include  not  only  mechanical 
displacements,  but  also  the  electrical  displacements  or  charge  flow  depending  on  how  the  system 
of  equations  was  developed. 

3.1.8.  System  reduction 

Often  it  is  desired  to  reduce  a  system  using  the  mode  shapes  rather  than  work  with  the 
full  finite  element  model.  This  drastically  reduces  the  number  of  degrees  of  freedom  and  allows 
for  much  faster  computations  once  the  eigen  system  is  solve  for.  This  is  often  done  in  structural 
mechanics  because  the  necessary  assumption  that  the  modal  matrix  (O)  uncouples  the  damping 
matrix  C  is  usually  valid  for  structural  damping.  That  is 

0>TCO  =  diag[u,  f22  •••  Hm]  (3-82) 

However,  this  assumption  cannot  be  made  for  most  electrical  systems  that  would  be  connected  to 
the  PZTs.  So,  in  general,  the  coupled  electrical-mechanical  system  of  Eq.  (3.78)  could  not  be 
reduced  in  its  entirety  using  a  few  of  its  mode  shapes. 
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The  system  can  be  partially  reduced  using  the  structural  mode  shapes  if  the  assumption  is 
made  that  the  mode  shapes  for  the  actual  coupled  system  can  be  expressed  in  terms  of  the  open 
circuit  mode  shapes.  This  assumption  is  generally  very  reasonable  since  changes  in  the  elastic 
stiffness  of  PZTs  due  to  open  circuiting  or  short  circuiting  cause  only  modest  shifts  in  natural 
frequency  and  very  small  alterations  in  the  mode  shapes.  The  structural  stiffness  added  by  the 
PZTs  is  generally  small  compared  to  the  stiffness  of  the  base  structure  and  the  difference  between 
open  circuit  and  short  circuit  stiffness  is  only  around  ten  to  twenty  percent.  This  combined  with 
the  fact  that  PZTs  usually  only  cover  a  fraction  of  the  structural  surface  means  that  the  differences 
in  mode  shapes  between  the  open  circuited  condition  and  the  actual  structure  will  be  generally 
small  and  localized. 

First,  the  coupled  system  is  reduced  so  that  the  open  circuit  eigenvalue  problem  can  be 
solved  for  the  desired  number  of  eigenvalues. 

Kuu<p  =  C02ocMa<p  (3.83) 

Then  using  Eq.  (3.82)  and 

®TK„„®  =  diag[i»il  -  ®L]  as4' 

4>TM„4>  =  Im  (3.85) 

the  coupled  system  of  Eq.  (3.78)  can  be  reduced  to 

X  °lffl  Jdiagkl  °lfrl  TdiagR]  ®Tk«,1{ r  1  =  foX]  (3.86) 

_  0  MqJ[qJ  [  0  CqJ\qeJ  Kqu<D  K;q  J\qJ  [  Fq  J 

where 

r  =  0-1ue  (3.87) 

Now  the  problem  has  been  reduced  to  a  small  system  composed  of  only  the  electrical  degrees  of 
freedom  and  the  chosen  number  of  mode  shapes. 
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3.2.  Results  and  Verification 

In  this  section  a  comparison  is  made  between  the  developed  model  and  known  solutions 
to  problems.  First  the  finite  element  model  is  validated  by  comparing  static  deflections  and 
natural  frequencies  with  published  results  for  both  isotropic  and  orthotropic  plates.  Next  the 
results  predicted  by  the  model  for  plates  with  passive  electrical  damping  circuits  are  compared 
with  experimental  results,  demonstrating  the  ability  of  the  developed  theory  to  model  smart 
structural  systems. 


Table  3.1. 

Finite  element  convergence. 


Finite  element  mesh  size 

Result 

4x4 

8x8 

16x16 

32x32 

64x64 

w 

0.04729 

0.04512 

0.04457 

0.04443 

0.04440 

(J\ 

0.31304 

0.29353 

0.28887 

0.28772 

0.28743 

3.2.1.  Finite  Element  Verification 

In  order  to  verify  correct  finite  element  modeling  and  implementation  of  the  higher  order 
laminate  theory  a  test  specimen  is  modeled  to  allow  comparison  with  published  results.  A 
rectangular  plate  simply  supported  on  all  four  sides  is  modeled  under  a  uniform  pressure.  The 
plate  has  edges  of  lengths  a  and  b,  and  has  thickness,  h.  The  displacements  and  stresses  have  been 
normalized  using  the  following  relation 


<7  i  =  a  i 


-  fa  b)h3E2 
w- w  — ,  4 

U  2  )q0a 

(3.88a) 

a  b  h\  h2  .  10 

— ,± - r  i=l,2 

o  9  n  9  o  ~  ~ 2 

(3.88b) 
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—  (  h\  h2 

a e  =  cr6  0,0, +—  - j  (3.88c) 

V  2J<loa 

S.=a4(£.0,o)—  (3.88d) 

U  )q0a 

a!=0'/o,-,o'|—  (3.88e) 

V  2  Jq0a 

The  first  check  that  was  performed  on  the  finite  element  model  was  to  check  that  the 
results  converged  with  an  increasing  number  of  nodes.  The  ratio  of  the  length  of  the  plate  edge,  a, 
to  the  thickness,  h,  was  taken  to  be  100  and  the  plate  was  assumed  to  be  square  (a=b).  The  plate 
is  assumed  to  be  isotropic  with  a  Poisson’s  ratio,  v,  of  0.3  and  an  elastic  modulus,  E,  of  160Gpa. 
Table  3.1  shows  the  out-of-plane  displacement  and  the  maximum  in-plane  stress  at  the  center  of 


Table  3.2. 

Stress  and  displacement  comparison  for  an  isotropic  simply  supported  plate. 


5 

10 

100 

CPT 


Source 

w 

±Gi 

±g2 

06 

-ct4 

-a5 

FEM 

0.5369 

0.2944 

0.2944 

0.2210 

0.4956 

0.4956 

Reddy 

0.0535 

0.2944 

0.2944 

0.2112 

0.4840 

0.4840 

FEM 

0.04679 

0.2889 

0.2889 

0.2031 

0.5031 

0.5031 

Reddy 

0.0467 

0.2890 

0.2890 

0.1990 

0.4890 

0.4890 

FEM 

0.04451 

0.2871 

0.2871 

0.1957 

0.5062 

0.5062 

Reddy 

0.0444 

0.2873 

0.2873 

0.1947 

0.4909 

0.4909 

FEM 

0.04449 

0.2871 

0.2871 

0.1956 

- 

- 

Reddy 

0.0444 

0.2873 

0.2873 

0.1946 

- 

- 

FEM 

0.1251 

0.6200 

0.2815 

0.3058 

0.6896 

0.5455 

Reddy 

0.1248 

0.6202 

0.2818 

0.2927 

0.6745 

0.5201 

FEM 

0.1144 

0.6121 

0.2783 

0.2863 

0.6955 

0.5512 

Reddy 

0.1142 

0.6125 

0.2789 

0.2809 

0.6794 

0.5230 

FEM 

0.1109 

0.6095 

0.2773 

0.2789 

0.6977 

0.5533 

Reddy 

0.1106 

0.6100 

0.2779 

0.2769 

0.6813 

0.5240 

FEM 

0.1109 

0.6094 

0.2773 

0.2788 

- 

- 

Reddy 

0.1106 

0.6100 

0.2779 

0.2769 

- 

- 
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the  plate  using  differing  numbers  of  elements.  The  results  show  good  convergence  with  an 
increasing  number  of  nodes  for  both  displacements  and  stresses,  and  the  results  using  a  modest 
number  of  elements  per  side  (~10-20)  are  seen  to  be  reasonably  accurate. 

Next  the  finite  element  implementation,  of  the  higher  order  plate  theory  is  verified  against 
the  numerical  results  presented  by  Reddy  (1984b).  The  plate  under  consideration  is  isotropic 
(v=0.3)  and  comparison  is  made  for  a  variety  of  plate  thicknesses.  The  finite  element  model  used 
19  by  19  elements  for  the  square  plate  and  13  by  25  elements  for  the  rectangular  plate.  The 
results  are  presented  in  Table  3.2. 

Next  a  comparison  is  made  using  an  orthotropic  material.  The  plate  had  an  edge  length 
of  6=0. 32m,  and  a  IMpa  uniform  pressure  was  applied.  For  the  case  of  b/a  =1  a  19  by  19 
element  mesh  was  used.  For  the  case  of  b/a=2  a  13  by  25  element  mesh  was  used.  The 
orthotropic  material  properties  used  were 

Ex=  143.6GPa  £2  =75.43GPa  G,2  =42.06GPa  (3  g9) 

G13  =  25.58GPa  G23  =42.68GPa  v12=0.44  v21=0.23 

The  results  for  the  finite  element  model  as  well  as  the  exact  elasticity  solution  from  Srinivas,  Rao 

and  Rao  (1970)  are  shown  in  Table  3.3. 

A  final  verification  is  made  using  the  computed  natural  frequencies  for  a  square  plate. 


Table  3.3. 

Stress  and  displacement  comparison  for  an  orthotropic  plate. 


b/a 

h/a 

w  (mm) 

0,  (Mpa) 

a5  (Mpa) 

Exact 

HOT 

CPT 

Exact 

HOT 

CPT 

Exact 

HOT 

CPT 

0.05 

1.044 

1.045 

1.029 

144.3 

144.6 

144.2 

10.85 

11.17 

1 

0.10 

0.1377 

0.1369 

0.1286 

36.01 

36.45 

36.06 

5.382 

5.562 

- 

0.14 

0.05350 

0.05282 

0.04686 

18.34 

18.79 

18.40 

3.805 

3.953 

0.05 

1.077 

1.074 

1.064 

262.7 

262.4 

262.0 

14.05 

14.32 

- 

2 

0.10 

0.1408 

0.1384 

0.1329 

65.97 

65.90 

65.49 

6.927 

7.138 

0.14 

0.05421 

0.05235 

0.04845 

33.86 

33.82 

33.41 

4.878 

5.082 

- 
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The  plate  is  assumed  to  be  made  of  an  isotropic  material  with  a  Poisson’s  ratio,  v,  of  0.3,  an 
elastic  modulus,  E,  of  160Gpa  and  a  density  of  2000kg/m3.  A  19  by  19  element  mesh  is  used  for 
modeling,  and  the  length  to  thickness  ratio  ( a/h )  considered  is  10.  The  results  are  presented  in 
Table  3.4  along  with  the  exact  values  from  Srinivas,  Rao  and  Rao  (1970).  The  natural 
frequencies  have  been  normalized  using 


0)  =  C0 


(3.90) 


The  results  show  good  correlation  with  the  exact  solution  values.  The  differences  seen  in 
the  higher  frequency  modes  are  a  result  of  the  coarseness  of  the  finite  element  mesh  and  a  finer 
mesh  size  would  provide  more  accurate  results. 


Table  3.4. 


Natural  frequencies  computed  using  finite  element  model  compared  with  exact  results. 


m 

n 

Exact 

HOT 

CPT 

i 

i 

0.0932 

0.09289 

0.09539 

2 

1 

0.226 

0.2214 

0.2354 

2 

2 

0.3421 

0.3386 

0.3710 

1 

3 

0.4171 

0.4140 

0.4617 

2 

3 

0.5239 

0.5166 

0.5903 

1 

4 

- 

0.6508 

0.7649 

3 

3 

0.6889 

0.6751 

0.7987 

2 

4 

0.7511 

0.7386 

0.8846 

3.2.2 .  Passive  Electrical  Damping 

To  study  the  validity  of  simultaneously  modeling  the  electrical  and  mechanical  aspects  of 
a  smart  structural  system  it  is  necessary  to  compare  the  model  results  with  experimental  data  for  a 
structural  system  with  a  well-defined  electrical  system.  One  such  smart  structural  system  is  that 
of  a  piezoelectric  patch  connected  to  a  passive  electrical  damping  circuit.  Passive  damping 
circuits  usually  consist  of  only  inductors,  resistors  and  sometimes  capacitors,  making  them  easy 
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Fig.  3.2.  Top  view  of  cantilevered  plate  with  driving  and  shunted  piezoceramic  pairs. 

to  model  with  linear  electrical  equations.  Also,  there  is  a  far  amount  of  literature  that  has  been 
published  on  passive  damping  circuits,  providing  a  good  deal  of  experimental  data  ideal  for 
model  validation. 

A  cantilevered  plate  with  a  shunt  damping  circuit  studied  by  Hagood  and  Von  Flotow 
(1991)  is  first  examined.  The  aluminum  plate  seen  in  Fig.  3.2  has  two  pairs  of  PZTs  mounted  on 
the  upper  and  lower  surfaces.  One  pair  is  actuated  in  opposite  directions  to  drive  the  vibration  of 
the  plate  and  the  other  is  connected  to  an  inductor  and  resistor  in  series.  The  material  properties 
of  the  plate  and  PZTs  are  listed  in  Table  3.5. 

The  cantilevered  plate  was  modeled  with  a  28x3  finite  element  mesh.  The  higher  order 
laminate  theory  results  in  7  degrees  of  freedom  per  node  and  812  degrees  of  freedom  total.  The 


Table  3.5. 

Material  properties  from  experiments  of  Hagood  and  Von  Flotow. 


Material  property 

Aluminum 

PZT 

E  (Gpa) 

73.0 

63.0 

V 

0.33 

0.3  a 

p  (kg/m3) 

2800  a 

7750  a 

d„  (m/V) 

n/a 

180xl0'12 

e3s  (F/m) 

n/a 

9.22xl0"9 

Assumed  values 
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charge  across  the  four  PZT  patches  was  described  by  four  additional  degrees  of  freedom.  The 
open  circuit  mass  matrix,  Muu,  and  stiffness  matrix,  Kuu,  are  used  to  solve  for  the  first  twelve 
open  circuit  modes  with  a  LAPACK  generalized  eigen  solver.  The  system  is  then  reduced  using 
Eq.  (3.86)  in  order  to  determine  the  frequency  response.  The  electrical  equations  for  the  each 
circuit  are  now  added  directly  to  the  resulting  set  of  equations.  A  damping  coefficient  of  £=0.01 
was  assumed  for  each  mode  to  create  the  damping  matrix, 

The  model  is  first  analyzed  using  an  untuned  shunt  circuit  containing  only  a  resistor  and 
no  inductor.  The  increase  in  the  damping  coefficient,  £,  for  the  first  bending  mode  is  then 
calculated  for  a  variety  of  resistance  values.  These  results  are  shown  in  Fig.  3.3  against  the  non- 
dimensional  resistance.  The  non-dimensional  resistance  is  the  ratio  of  the  actual  resistance  to  the 
optimum  value  of  28.68kD  predicted  by  Hagood  and  Flotow  (1991).  The  results  show  good 
correlation  with  the  published  experimental  values.  Since  the  structural  damping  of  the  beam 
with  the  PZTs  shorted  out  was  only  0.16%,  the  untuned  shunt  damping  circuit  can  be  seen  to 
roughly  quadruple  the  damping  of  the  structural  system. 


Fig.  3.3.  Added  damping  to  first  mode  of  the  cantilevered  plate  for  untuned  resistance  shunting. 


Fig.  3.4.  Analytical  strain  for  first  mode  actuation  of  cantilevered  plate  with  shunted 
piezoceramic  pair  open,  shunted  and  optimally  tuned  and  damped. 


The  plate  is  then  modeled  with  a  tuned  shunt  damping  circuit.  The  inductor  used  has  the 
optimum  value  of  148.2  H  predicted  by  the  literature  (Hagood  and  Flotow,  1991).  It  is  connected 
in  series  with  a  resistor  having  a  varying  amount  of  resistance.  Figure  3.4  shows  the  first  bending 
mode  frequency  response  for  the  plate  with  three  conditions  for  the  shunted  PZTs;  open  circuited, 
shorted  out,  and  optimally  damped  using  the  values  of  Hagood  and  Flotow  (1991).  The  results 
clearly  show  the  reduction  in  natural  frequency  caused  by  shorting  out  the  PZT.  Also  Fig.  3.4 
shows  how  efficient  the  tuned  damping  circuit  can  be  at  damping  out  a  particular  mode  of 
vibration.  The  frequency  response  for  the  first  mode  is  also  shown  in  Fig.  3.5  for  varying 
resistance.  When  the  resistance  is  infinitely  large,  the  system  behaves  as  in  open  case.  As 
resistance  decreases,  the  damping  increases,  flattening  the  system  response.  As  resistance 
decreases  below  the  optimum  value  of  6640Q  the  system  begins  to  develop  two  separate  peaks. 


t 


47 

These  results  show  fair  correlation  with  the  experimental  data  of  Hagood  and  Flotow  (1991) 
presented  in  Fig.  3 .6;  however,  there  are  still  noticeable  differences.  This  is  most  likely  the  result 
of  inaccuracy  in  the  predicted  natural  frequency  of  the  current  theory.  The  finite  element  model 
predicts  an  open  circuit  natural  frequency  of  34.136  Hz  as  opposed  to  the  measured  natural 
frequency  of  33.77  Hz.  This  is  most  likely  caused  by  a  combination  of  the  finite  element 
discretization  and  difference  between  the  assumed  material  properties  and  the  actual  properties. 
When  the  finite  element  model  is  modified  using  lower  stiffnesses  to  closely  match  the 
experimental  value  of  natural  frequency,  the  frequency  responses  are  almost  identical  to  the 
experimental  results  of  Hagood  and  Flotow  (1991),  as  shown  in  Fig.  3.7. 


Fig.  3.5.  Analytical  strain  for  first  mode  actuation  of  cantilevered  plate  with  tuned  shunt  circuit 


and  varying  resistance. 
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Fig.  3.6.  Hagood’s  experimental  results  for  first  mode  actuation  of  cantilevered  plate  with  tuned 

shunt  circuit  and  varying  resistance. 


Fig.  3.7.  Strain  for  first  mode  actuation  of  the  modified  cantilevered  plate  with  tuned  shunt  circuit 


and  varying  resistance. 


49 


* 


Fig.  3.8.  Frequency  response  for  aluminum  plate  with  mode  2  damped. 


Next  the  multi-mode  damping  technique  of  Wu  (1998)  is  used  on  the  same  cantilevered 
aluminum  plate.  When  a  single  mode  is  damped,  an  inductor  and  a  resistor  in  parallel  are 
connected  to  the  shunted  PZTs.  This  leads  to  electrical  equations  of  the  form 


L 

-L 


-L 

L 


q  + 


o  o 
0  R 


q+ 


o  o 
0  0 


q  = 


(3.91) 


where  L  and  R  are  the  inductance  and  resistance  respectively.  These  are  then  added  to  the 
structural  finite  element  equations.  The  values  of  the  inductance  and  resistance,  determined  using 
the  method  proposed  by  Wu  (1998),  are  199.4H/113.7k£2,  8.888H/27.66kn,and  1.188H/16.38kf) 
for  bending  modes  1,  2,  and  3  respectively.  The  tip  deflection  as  a  function  of  frequency  for  the 
plate  with  mode  2  damped  is  plotted  in  Fig.  3.8  with  that  of  the  shorted  plate.  Next  modes  1  and 
2  are  damped  with  the  circuit  shown  in  Fig.  3.9.  The  values  of  L2  and  C2  are  chosen  such  that 
L'2C2  =  1 Jco2x  and  are  assumed  to  be  22.5 16H  and  l.OOpF  respectively.  The  values  of  R2"  and 
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I2”,  calculated  using  the  multi-mode  theory  of  Wu  (1998),  are  36.54kQ  and  9.862H  respectively. 
This  leads  to  the  electrical  equations  given  by 

r-  ,  r  „  nr.  [0000  0 

0  0  0  0  0  R,  -/?,  0  0  !  -i 

0  0  0  0  0  -fl,  R,  +  Rl  Rt  0  -R’  °  ~C2  °  ~C2 

0  0  I,  0  0  q+  -Rt  R}  0  0  q+  0  0  0  0  0 

o  o  o  l;  o  o  o  ooo  0  ~c  0  ~c  0 

_0  0  0  0  L'J  [  o  -K  0  0  K  \  [ooo  0  0 

The  resulting  frequency  response  is  shown  in  Fig.  3.10.  This  technique  can  be  seen  to  result  in  a 
decrease  in  maximum  tip  deflection  of  over  one  order  of  magnitude.  This  versatile  technique  can 
also  be  extended  to  as  many  modes  as  one  wishes  to  damp  out. 

Finally  a  comparison  is  made  between  the  open  circuited  and  short  circuited  mode 
shapes.  As  discussed  above,  the  system  can  be  reduced  using  its  open  circuited  eigenvectors. 
However,  it  is  desirable  to  estimate  the  nature  of  any  resulting  errors  when  this  technique  is 
applied  to  an  electrical  system  where  the  short  circuit  mode  shapes  are  more  appropriate.  Since 
the  open  and  shorted  eigenvectors  are  both  normalized  to  the  same  mass  matrix,  the  difference 
between  the  two  can  be  directly  calculated.  Then  to  compare  the  deviation,  the  norm  of  the 
difference  can  be  compared  to  the  norm  of  the  original  eigenvector.  The  ratio  of  the  2-norms  of 

tfl  qi 

PZT 


Fig.  3.9.  Two  mode  damping  circuit. 
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Fig.  3.10.  Frequency  response  for  aluminum  plate  with  modes  1  and  2  damped. 


the  difference  to  the  open  circuit  mode  shape  is  found  to  be  1.79x10 s,  4.53x10 5,  and  2.73x10 
for  the  first  three  bending  modes  of  the  cantilevered  aluminum  plate  use  above.  These  very  small 
differences  are  negligible  when  compared  to  the  changes  in  natural  frequency,  which  are  on  the 
order  of  one  percent.  Thus,  the  technique  for  using  open  circuit  mode  shapes  for  system 
reduction  will  likely  yield  good  results  in  most  practical  problems. 

3.3.  Transient  Analysis 

Although  the  dynamic  results  shown  thus  far  are  useful  for  the  study  of  smart  structures, 
the  ability  to  model  the  transient  response  of  such  structures  is  the  most  important  and  most 
difficult  goal  to  achieve.  Active  control  systems  are  often  used  to  damp  out  transient  vibrations 
using  feedback  control  circuits,  thus  an  accurate  analysis  of  the  control  system  requires  an 
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accurate  estimation  of  the  dynamic  response  of  the  system.  Also,  damage  detection  schemes 
often  use  changes  in  the  dynamic  response  of  the  structure  to  detect  and  localize  damage.  Thus, 
an  accurate  modeling  technique  is  an  invaluable  tool  in  the  design  such  systems.  Having 
demonstrated  the  basic  validity  of  the  developed  theory,  it  is  desired  to  study  its  ability  to 
accurately  model  the  transient  response  of  adaptive  structural  systems  and  compare  the 
mechanical  arid  electrical  response  with  experimental  results  and  other  modeling  techniques. 

3.3.1.  Newton-Raphson  Method 

The  nonlinear  transient  analysis  is  conducted  using  Newmark-beta  method  with  Newton- 
Raphson  (N-R)  iteration  (Argyris  and  Mlejnek,  1991;  Cook,  et  al.,  1989;  Bathe,  1996).  The  above 
system  of  dynamic  equilibrium,  Eq.  (3.81),  can  be  rewritten  in  the  implicit  form  at  time  (t  +  AO  as 

follows 

ma(i  +  A*)*  +  c (t  +  A t)v(t  +  A t)k  +  k(/  +  At)Au*  =  f(£  +  A t)  -  R(t  +  At)  (3.93) 
where,  m  is  the  time  independent  augmented  mass  matrix,  c(t+At)  is  the  augmented  damping 
matrix  at  time  (t+At)  and  k(t+At)  is  the  augmented  stiffness  matrix  at  time  (t+At).  The  vector, 
f(t+At),  contains  the  externally  applied  nodal  point  loads,  which  is  the  right  hand  force  vector  in 
Eq(3.81)  and  the  vector,  R(t+At)(k'1),  contains  nodal  point  forces  equivalent  to  elemental  stresses 
at  time  (t+At)  in  (k-l)th  iteration.  The  mechanical  and  electric  displacements  at  the  nodes  are 
contained  in  the  vector,  u. 

The  incremental  nodal  displacements  at  the  kth  iteration  of  a  particular  time-step  (t+At) 
are  defined  for  a  given  time  step,  At,  as  follows 


Au*  =  u(/  +  A t)k  -  u (t  +  A /)*“' 
Using  the  trapezoidal  rule,  the  following  assumptions  are  made 


(3.94) 
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u  (t  +  At)  =  u(0  +  —  {v(0  +  v(f  +  AO} 


(3.95a) 


y(t  +  At)  =  \(t)  +  —  (a(0  +  a(?  +  A0} 


(3.95b) 


This  leads  to  the  following  forms 


y(t  +  At)  =  {u(t  +  At)  -  o(0}--  v(0 


(3.96a) 


4  4 

a(t  +  At)  =  {u  (t  +  At)  -  u(0}  - — v(0  “  a(0 


(3.96b) 


For  the  k*  iteration  of  a  particular  time-step  (t+At)  substitution  of  Eq.  (3.94)  into  Eqs.  (3.96a, b) 


results  in  the  following 


y(t  +  At)k  =  -^{u(t  +  At)k  ' +Au* -u(0}-— v(0 


(3.97a) 


a(t  +  At)k  =  — ^—r  (u(^  +  A?)*-1  +Au*  -ii(0}— T7v(0_a(0 
V  7  (A/)2  A* 


(3.97b) 


Substitution  of  Eqs.  (3.97a,b)  into  Eq.  (3.93)  yields  the  iterative  form 


kAu*  =  f(t  +  At)-R(t  +  At)k~'  -maf"1  -c(t  +  At)a\ 


(3.98) 


where 


a}-1  =[— ^y{u(t  + AO*'1  -u(0}“~v(0-a(0] 


(3.99a) 


a*'1  =  [—  {u(/  +  AO*'1  -  u(0}  -  v(0] 


(3.99b) 


A  O 

k  = - -m  +  —  c(t  +  At)  +  k(t  +  At) 

(At)2  At 


(3.99c) 


R  (t  +  AO*'1  =  k(/  +  At)u(t  +  AO* 


(3.99d) 


This  results  in  a  time  integration  method  that  can  be  iterated  at  each  time  step  to  provide  accurate 


prediction  of  the  transient  response  of  the  system. 


t 
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1.83  cm 


5.1  cm 


3.3.2.  Results 

The  developed  model  is  used  to  calculate  the  response  of  a  composite  plate  with  surface 
bonded  actuators  subjected  to  impulse  loading.  The  objective  is  to  demonstrate  the  nature  and 
magnitude  of  errors  that  exist  when  simpler  approaches  are  used.  First,  a  comparison  is  made 
between  the  response  predicted  by  the  refined  higher  order  laminate  theory  and  the  classical  plate 
theory.  Then  the  difference  between  the  coupled  piezoelectric-mechanical  model  and  the 
uncoupled  model  is  examined. 

The  plate  under  consideration  is  clamped  at  one  end  with  a  cantilevered  section  31.1cm 
long  by  5.1cm  wide.  The  plate  is  a  graphite-epoxy  laminate  with  sixteen  plies  of  0.137mm  ply 
thickness.  A  variety  of  ply  stacking  sequences  are  considered,  including  cross-ply  as  well  as 
balanced  and  unbalanced  angle-ply  lay-ups.  A  piezoelectric  patch  is  bonded  to  the  upper  surface 
of  the  plate  as  shown  in  Figure  3.11.  The  patch  is  made  of  PZT-5H  with  a  thickness  of  0.25mm. 
The  material  properties  for  the  graphite-epoxy  and  PZT  are  listed  in  Table  3.6. 
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Table  3.6. 

List  of  material  properties  used  for  transient  analysis. 


Property 

Graphite-epoxy 

PZT-5H 

E,  (GPa) 

372.0 

60.6 

E2  (GPa) 

4.12 

60.6 

V 12 

0.275 

0.29 

V23 

0.42 

0.48 

G,2  (GPa) 

3.99 

23.49 

G23  (GPa) 

3.6 

22.99 

p  (kg/m3) 

1788.5 

7500 

d3|  and  d32  (nm/V) 

- 

-0.274 

d24  and  d)5  (nm/V) 

- 

0.741 

XT33  (nF/m) 

- 

30.1 

Fig.  3.12.  Tip  displacement  for  [0°,90o]4s  laminate  (L/h=142)  under  lps  impulse  loading. 


3.3.2. 1  Higher  order  versus  classical  plate  theory 


First,  a  comparison  is  made  between  the  refined  higher  order  laminate  theory  and  the 
classical  plate  theory.  A  laminate  stacking  sequence  of  [0°,90°]4S  is  considered  first,  with  the 
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PZT  left  open  circuited.  This  laminate  has  a  ratio  of  plate  length  to  thickness  of  142,  making  this 
a  relatively  thin  plate.  The  plate  is  modeled  with  the  refined  higher  order  laminate  theory,  which 
includes  the  effects  of  transverse  shear,  as  well  as  with  the  classical  plate  theory,  which  neglects 
transverse  shear.  All  other  aspects  of  the  modeling  are  identical  in  both  cases.  The  plate  is 
subjected  to  a  5-Newton,  lps  impulse  point  load  at  the  tip.  The  analysis  is  performed  for  both 
cases  with  a  lps  time  step.  Figure  3.12  shows  the  resulting  tip  displacement  obtained  using  the 
two  models.  The  difference  between  the  higher  order  theory  and  the  classical  theory  is  modest, 
but  noticeable  during  the  short  time  interval  analyzed,  even  for  this  thin  laminate.  By  including 
transverse  shear  stress,  the  refined  higher  order  theory  results  in  a  model  with  lower  natural 
frequencies  than  those  predicted  by  the  classical  plate  theory.  This  reduction  in  natural 
frequencies  is  due  to  the  fact  that  during  vibration  the  plate  undergoes  not  only  pure  bending,  but 
also  deformation  in  the  form  of  out-of-plane  shear.  The  bending  modes  are  all  composed  of  both 


Fig.  3.13.  Sensor  output  for  [0°,90o]4s  laminate  (L/h=142)  under  lps  impulse  loading. 
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forms  of  deformation  and,  thus,  models  which  do  not  include  both  are  too  stiff.  The  effects  are 
not  only  more  significant  for  thicker  laminates,  but  also  for  the  higher  frequency  modes,  and  this 
difference  in  frequencies  can  be  seen  in  Figure  3.12.  Figure  3.13  shows  the  voltage  output  from 
the  PZT  sensor  during  the  impulse  loading.  A  very  significant  difference  in  voltage  output  is 
observed  between  the  two  plate  theories.  Although  the  shape  of  the  response  is  very  similar,  the 
magnitude  of  the  signal  differs  significantly.  This  is  because  at  these  high  frequencies  the  output 
of  the  coupled  piezoelectric-mechanical  theory  is  very  sensitive  to  changes  in  the  strain  field  of 
the  vibrational  waves  than  mode  across  the  sensor. 

Next,  the  same  plate  is  analyzed  under  a  1 -Newton,  5 (is  impulse  point  load  at  the  tip. 

The  time  step  for  the  transient  analysis  is  5 (is,  so  that  a  longer  time  segment  could  be  examined. 
The  tip  displacement  for  the  higher  order  and  classical  plate  theories  is  shown  in  Figure  3.14  and 
the  sensor  output  is  shown  in  Figure  3.15.  When  a  longer  time  step  is  used  the  differences 
between  the  two  theories  are  much  less  significant  for  the  thin  plate,  because  the  natural 


Fig.  3.14.  Tip  displacement  for  [0°,90o]4s  laminate  (L/h=142)  under  5ps  impulse  loading. 
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frequencies  associated  with  the  lower  order  modes  predicted  by  the  two  plate  theories  are  quite 
similar.  Note  that  in  this  case  the  charge  flow  from  the  PZT  is  measured  for  sensor  output.  It  can 
be  seen  that  the  differences  in  the  higher  frequencies  affect  sensor  output  to  a  greater  degree  than 
tip  displacement. 

An  increase  in  plate  thickness  increases  the  difference  in  the  response  predicted  by  the 
higher  order  theory  and  the  classical  plate  theory.  This  is  shown  in  Figure  3.16,  which  depicts  the 
tip  response  for  a  [0°5,90o5]4s  laminate.  The  ratio  of  plate  length  to  thickness  is  28.4  for  this  case, 
making  this  a  moderately  thick  plate.  The  impulse  length  and  time  step  used  in  this  case  are  both 
lps.  By  increasing  plate  thickness  differences  are  observed  even  in  the  low  frequency  modes 
predicted  by  the  two  plate  theories.  This  results  in  sizable  differences  in  plate  deflections. 


Fig.  3.15.  Sensor  output  for  [0°,90o]4s  laminate  (L/h=142)  under  5ps  impulse  loading. 
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33.2.2  Coupled  versus  uncoupled  piezoelectric  theory 

Next  a  comparison  is  made  between  the  coupled  and  the  uncoupled  piezoelectric  models. 
The  plate  with  [0°,90o]45  stacking  sequence  is  considered.  The  plate  is  subjected  to  a  1 -Newton, 
5|is  impulse  point  load  at  the  tip  and  a  5ps  time  step  is  used  for  the  transient  analysis.  The  PZT  is 
open  circuited  and  voltage  is  measured  for  sensor  output.  The  response  is  calculated  using  both 
the  coupled  theory  presented  in  this  work  as  well  as  the  traditional  uncoupled  approach.  The 
resulting  tip  displacement  is  shown  in  Figure  3.17  for  both  approaches.  Only  slight  differences 
between  the  two  models  are  observed  in  this  case.  This  is  due  to  the  relatively  small  size  of  the 
PZT  patch  and  its  limited  contribution  to  the  overall  plate  stiffness.  This  difference  is  expected  to 
be  larger  for  thicker  PZT  patches. 
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Figure  3.18  shows  the  charge  flow  of  the  sensor  for  both  cases.  Here  it  can  be  seen  that 
the  coupled  theory  predicts  dramatically  different  results  from  the  uncoupled  approach,  even 
though  the  displacements  are  shown  to  be  similar.  The  reason  for  the  large  differences  m  the 
electrical  response  is  due  to  a  combination  of  two  factors.  The  first  is  the  fact  that  the  uncoupled 
theory  assumes  that  the  electric  field  is  constant  over  the  entire  area  of  the  PZT.  During  impulse 
loading  high  frequency  bending  waves  travel  across  the  length  of  the  plate,  as  shown  m  Figure 
3.19.  As  these  high  frequency  waves  move  across  the  PZT,  the  piezoelectric  material  is  subjected 
to  areas  of  local  compression  and  tension.  As  a  result,  the  electric  displacement  is  positive  in 
some  local  areas  and  negative  in  others  as  shown  in  Figure  3.20.  The  net  charge  output  is  thereby 
reduced  since  charge  merely  flows  from  one  region  of  the  patch  to  another.  When  the  PZT  is 
open  circuited,  the  charge  flow  within  the  PZT  still  occurs  making  the  voltage  output  very 
sensitive  to  local  strain.  As  shown  in  Figure  3.21,  the  strains  at  the  two  ends  of  the  PZT  are  very 
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dissimilar  during  the  impulse  loading.  Thus,  there  exists  a  continuously  changing  strain  gradient 
across  the  PZT  patch,  which  leads  to  the  difference  in  results  predicted  by  the  coupled  and 
uncoupled  models.  Figure  3.22  presents  the  output  if  the  voltage  of  the  PZT  is  measured,  as 


Fig.  3.18.  Sensor  charge  flow  for  [0°,90o]4s  laminate  under  impulse  loading. 


Fig.  3.19.  Plate  displacement  at  t=73.95  ps  during  impulse  tip  loading. 
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opposed  to  the  charge  flow.  Again  the  coupled  theory  predictions  greatly  differ  from  the 
uncoupled  approach  and  the  high  frequency  signal  can  be  seen  to  be  much  more  prominent  when 
voltage  is  measured. 


Fig.  3.20.  Electric  displacement  at  t=60.45  ps  during  impulse  loading. 
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The  second  factor  that  creates  differing  results  between  the  coupled  and  uncoupled 
approaches  is  the  existence  of  electrical  inductance,  both  in  the  electrical  circuit  and  within  the 
piezoelectric  material  itself.  This  inductance  acts  likes  mass  in  a  mechanical  system  and  prevents 
instantaneous  changes  in  the  charge  flow  of  the  electrical  system.  The  magnitude  of  the  electrical 
signal  created  by  high  frequency  vibrations  in  the  piezoelectric  material  are  thus  limited  in 
magnitude,  just  as  mass  limits  amplitude  of  vibrations  induced  by  forces  in  a  mechanical  system. 
The  uncoupled  approach  is  incapable  of  including  inductance  effects,  thus  the  high  frequency 
component  of  the  electrical  sensor  output  will  always  be  overpredicted  when  compared  to  a 
physical  system.  The  coupled  model  has  the  ability  to  include  inductance  that  exists  in  the 
physical  system.  A  typical  lm  long  coax  cable  has  inductance  on  the  order  of  0.3p.H  and  this  can 
be  included  by  use  of  Eq.  (3.80).  The  piezoelectric  device  itself  also  has  some  inductance  which 
affects  not  only  the  net  output  of  the  device,  but  also  the  charge  flow  from  one  area  of  the 


Figure  3.22.  Sensor  output  for  [0°,90o]4s  laminate  under  impulse  loading. 


actuator  to  another.  This  inductance  can  be  included  in  Eq.  (3.75)  by  addition  of  positive  values 
along  the  diagonal  of  the  section  of  the  mass  matrix  corresponding  to  the  piezoelectric  electric 
displacement.  The  exact  value  for  the  internal  inductance  is  very  difficult  to  quantify  and 
comparison  with  experimental  data  seems  to  be  the  best  method  for  obtaining  these  values. 

3.3.3.  Experimental  comparison 

To  validate  the  model  and  verify  the  conclusions  made  from  comparison  of  the  coupled 
theory  with  the  uncoupled  approach,  a  set  of  experiments  was  performed.  The  objective  of  these 
experiments  was  to  provide  a  set  of  transient  sensor  outputs  for  a  well-defined  model  and 
compare  these  directly  with  the  results  predicted  by  the  developed  model. 

The  test  specimen  used  was  a  cantilevered  plate  with  two  piezoelectric  patches  bonded  to 
the  surface.  The  plate  was  made  from  3.18mm  thick  Aluminum  2024-T3  and  the  plate  geometry 


Fig.  3.23.  Experimental  setup  and  plate  geometry  for  transient  analysis. 
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is  that  shown  in  Figure  3.23.  The  small  piezoelectric  patch  at  the  root  was  a  custom  ACX 
actuator  with  wafer  dimensions  of  1.27cm  length,  0.635cm  width  and  0.25mm  thickness.  The 
larger  piezoelectric  patch  was  an  ACX  QP-10N  actuator  with  wafer  dimensions  of  4.597cm 
length,  2.057cm  width  and  0.25mm  thickness.  Both  devices  were  made  of  PZT-5A  with  a 
polyamide  coating.  The  material  properties  for  the  aluminum  and  the  piezoelectric  devices  are 
shown  in  Table  3.7. 

The  small  patch  at  the  root  was  used  as  an  actuator  to  induce  vibration  in  the  plate,  while 
the  larger  patch  was  used  as  a  sensor.  To  simulate  an  impulse  loading  in  the  plate,  a  single  cycle 
of  high  frequency  voltage  was  applied  to  the  actuator.  The  single  cycle  voltage  was  generated  by 
a  function  generator  set  operate  in  burst  mode  and  create  a  3.75V  signal  once  every  two  seconds. 
The  input  voltage  was  then  amplified  to  a  75V-peak  signal  that  was  sent  to  the  actuator.  This 
signal  is  illustrated  in  Fig.  3.24. 

The  sensor  output  was  measured  using  two  different  methods.  First  the  sensor  was 
connected  to  a  charge  amplifier  that  converted  to  the  charge  flow  of  the  device  into  a  voltage 
signal  that  was  measured  and  analyzed  using  an  oscilloscope.  This  method  allowed  measurement 
of  the  charge  output  of  the  piezoelectric  device  for  comparison  with  the  developed  model.  Next, 


Table  3.7. 

Material  properties  for  experimental  comparison. 

Material  property  Aluminum 


E  (Gpa) 

68.5 

60.6 

V 

0.326 

0.29 

G  (Gpa) 

25.83 

23.49 

p  (kg/m3) 

2784 

7500 

d3i,  d32  (m/V) 

n/a 

-274x1  O' 

d24,  du  (m/V) 

n/a 

741x10'’ 

e3s  (F/m) 

n/a 

3.01x10' 
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Fig.  3.24.  Input  voltage  used  to  simulate  impulse  loading. 

the  sensor  was  connected  directly  to  the  oscilloscope  and  the  voltage  output  was  measured. 

Sensor  output  was  measured  for  a  range  of  input  frequencies.  The  frequency  of  the  single 
cycle  voltage  input  allowed  control  of  the  duration  of  the  impulse.  Higher  frequency  inputs  acted 
like  shorter  duration  impulse  loads  and  excited  higher  frequency  responses  from  the  system.  This 
method  proved  to  be  much  more  controlled  than  mechanical  impulse  loads  and  the  output  signals 
for  the  impulses  were  very  repeatable.  To  minimize  any  effects  of  noise  and  deviation  between 
individual  impulses,  the  oscilloscope  was  set  to  provide  an  averaged  output  using  sixteen  impulse 

signals. 

The  charge  output  for  10kHz,  25kHz  and  50kHz  input  frequencies  is  shown  in  Fig.  3.25. 
The  important  aspect  to  be  noted  from  these  results  is  that  the  overall  shape  of  the  output  is  the 
same  for  all  three  frequencies.  The  effect  of  electrical  inductance  discussed  above  can  be  clearly 
seen,  in  that  the  magnitude  of  the  very  high  frequency  component  of  the  output  is  much  smaller 
than  the  lower  frequency  component,  which  governs  the  overall  shape  of  the  output.  Although 
the  higher  frequency  input  voltages  provide  more  high  frequency  excitation  of  the  system,  the 
magnitude  of  the  electrical  signal  created  by  this  component  is  limited  due  to  the  electrical 
inductance  and  the  fact  that  the  high  frequency  vibrations  cause  less  net  charge  flow. 
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Figure  3.26  shows  the  voltage  output  for  1kHz,  10kHz  and  25kHz  input  frequencies.  It  is 
clear  that  the  voltage  is  more  sensitive  to  high  frequency  vibration  and  less  affected  by  the 
electrical  inductance  than  the  case  when  sensor  charge  flow  is  used.  This  is  reasonable  since 
although  the  oscilloscope  is  not  a  true  open  circuit  it  has  a  1MQ  internal  resistance  that  limits  the 
electrical  system  to  very  small  net  charge  flows.  Less  charge  flow  in  turn  implies  less  effect  from 
inductance,  but  charge  flow  within  the  piezoelectric  device  during  high  frequency  vibration  is  still 
affected. 

The  developed  analysis  is  now  used  to  model  the  system  described  above.  A  finite 
element  mesh  with  85  by  9  elements  was  used  for  the  plate  resulting  in  860  nodes.  The  system 
had  6020  mechanical  degrees  of  freedom,  8  nodal  electrical  displacements  for  the  small  actuator 
and  48  nodal  electrical  displacements  for  the  sensor.  The  material  properties  listed  in  Table  3.7 
are  used  to  model  the  system.  The  sensor  charge  flow  predicted  by  the  model  for  10kHz  and 
25kHz  input  frequencies  are  shown  in  Figs.  3.27  and  3.28.  The  output  shows  good  correlation 
between  the  coupled  model  and  the  experimental  results  from  Fig.  3.25.  The  overall  shape  of  the 
signals  are  very  similar  and  similarities  between  the  two  frequencies  are  replicated  by  the  model. 
This  result  is  in  contrast  to  the  predicted  response  using  the  uncoupled  approach  seen  in  Fig.  3.29. 
The  lack  of  internal  inductance  causes  the  high  frequency  components  of  the  signal  to  dominate 
the  sensor  output.  The  results  from  the  coupled  approach  and  the  experimental  results  are  not 
identical  and  a  more  valuable  comparison  can  be  made  by  examining  the  response  in  the 
frequency  domain.  Figure  3.30  shows  the  frequency  content  of  the  response  for  the  coupled  and 
uncoupled  approaches.  Both  models  correlate  well  at  the  lower  frequencies,  but  the  uncoupled 
model  has  a  large  amount  of  high  frequency  content  as  discussed  above.  At  higher  frequencies, 
significant  differences  in  response  between  coupled  and  uncoupled  models  is  observed.  The 
results  from  the  coupled  model  also  deviates  from  the  experimental  results.  This  is  a  result  of  the 
limitations  of  the  finite  element  method  and  the  mesh  size  used  in  this  model.  The  variations 
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observed  are  not  unreasonable  considering  the  difficulties  associated  with  accurate  modeling  of 
high  frequency  structural  dynamics.  Figure  3.31  shows  a  comparison  between  the  classical  and 
higher  order  plate  models.  Both  results  were  computed  using  the  coupled  theory  and  show 
similar  response  curves.  The  primary  difference  between  these  two  models  appears  in  the  exact 
locations  of  the  response  peaks,  since  the  classical  plate  theory  specify  lower  natural  frequencies 
compared  to  the  higher  order  theory.  Once  again,  deviations  in  response  between  experiments 
and  models  are  observed  at  higher  frequencies  due  to  the  reasons  stated  above.  However,  the 
relative  merit  of  the  two  theories  is  difficult  to  establish  from  these  results. 

Next  sensor  voltage  was  calculated  using  the  developed  model.  The  sensor  voltage 
predicted  by  the  model  for  1kHz  and  10kHz  input  frequencies  are  shown  in  Figs.  3.32  and  3.33. 
The  output  is  similar  to  the  experimental  data,  but  is  less  accurate  than  the  predicted  charge  flow. 
The  reason  for  this  difference  is  the  difficulty  in  accurately  characterizing  the  electrical  system. 
The  oscilloscope  used  to  measure  the  voltage  experimentally  is  not  an  ideal  voltage-measuring 
device.  The  actual  circuit  measures  the  voltage  across  a  large  resistor  and  the  circuit  has  a  small 
amount  of  charge  flow,  contrary  to  the  zero  net  charge  flow  assumed  by  the  piezoelectric 
boundary  conditions.  Also,  the  charge  flow  within  the  PZT  becomes  extremely  important  in 
determining  the  output  voltage.  A  more  detailed  characterization  of  the  internal  inductance  and 
resistance  of  the  piezoelectric  device  and  voltage  measuring  circuit  must  be  used  to  accurately 
model  the  voltage  response  of  an  adaptive  structure.  The  sensor  charge  flow  is  less  sensitive  to 
these  effects,  thus  making  it  easier  to  obtain  an  accurate  model. 


Sensor  Output  (Coul) 


Fig.  3.3 1 .  Comparison  of  classical  and  higher  order  plate  theories  for  25kHz  input  frequency. 


Fig.  3.32.  Predicted  response  for  1kHz  input  frequency. 


4.  Delamination  Modeling  in  Composite  Laminates 


Delamination  occurs  when  a  region  of  composite  material  becomes  debonded  from  the 
layers  below  it.  Delamination  is  known  to  affect  not  only  the  compressive  strength  of  composite 
laminates,  but  also  the  dynamic  characteristics  of  a  composite  structure.  Since  adaptive  structural 
systems  using  piezoelectric  patches  are  often  used  to  detect  and  localize  damage,  such  as 
delamination,  the  developed  theory  is  extended  to  include  the  modeling  of  delamination.  The 
primary  goal  of  this  investigation  is  to  study  the  effects  of  delamination  on  the  system  response  as 
opposed  to  predicting  delamination  growth.  The  response  that  is  examined  is  not  just  the 
mechanical  response,  but  also  the  electrical  response  that  includes  the  sensor  signal.  Accurate 
prediction  of  the  sensor  signal  will  allow  future  damage  detection  studies  to  be  more  reliable. 

When  delamination  occurs,  a  thinner  sublaminate  is  formed  above  and  below  the  region 
of  delamination.  The  overall  method  for  delamination  modeling  used  in  this  work  is  to  model 
each  individual  sublaminate  as  a  separate  plate  and  then  join  the  sublaminates  together  to 
represent  the  entire  plate.  During  plate  vibration  the  possibility  exists  for  the  individual 
sublaminates  to  contact  or  impact  each  other,  so  a  method  for  analyzing  transient  response  that 
includes  this  contact  is  implemented.  During  transient  analysis  this  contact/impact  can  cause 
numerical  instabilities,  thus  a  discontinuous  time  integration  method  is  used  for  modeling  the 
transient  response  of  the  structure.  The  resulting  model  is  then  used  to  study  the  effect  of 
delamination  of  the  response  of  adaptive  composite  structures. 

4.1.  Modeling 

To  include  the  presence  of  delamination  the  laminate  theory  developed  in  the  previous 
chapter  must  be  modified  to  model  the  multiple  sublaminates  in  the  region  of  delamination.  This 
requires  modification  of  the  higher  order  displacement  field  and  also  selection  of  a  method  for 
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joining  the  sublaminates  together  at  their  boundaries  with  the  undamaged  host  structure.  In  this 
section  these  issues  are  addressed  and  the  model  including  delamination  is  developed.  It  should 
be  noted  that  delamination  growth  is  not  considered  in  this  work. 


4.1.1.  Higher  Order  Theory  with  Delamination 

As  discussed  in  the  previous  chapter,  the  refined  higher  order  laminate  theory  provides  a 
convenient  and  accurate  method  for  modeling  laminates  and  including  the  effect  of  transverse 
shear.  Transverse  shear  becomes  very  important  in  moderately  thick  plates,  and  as  a  result, 
classical  plate  theory  under  predicts  deformation  due  to  out-of-plane  loads  and  over  predicts 
natural  frequencies  for  the  bending  of  plates.  Delamination  creates  smaller  sublaminates  within 
the  host  structure.  Although  these  sublaminates  are  thinner  than  the  undamaged  laminate,  they 
are  often  relative  small  in  size  and  it  is  possible  that  the  sublaminates  will  have  a  higher  relative 
thickness  that  the  undamaged  laminate.  Thus,  transverse  shear  will  also  play  an  important  role 
for  accurate  modeling  of  local  vibration  modes  and  the  higher  order  laminate  theory  is  still  an 
effective  tool  for  addressing  transverse  shear  effects. 


Fig.  4.1.  Relationship  between  the  finite  element  mesh  of  the  sublaminates  and  the 
undelaminated  regions  of  the  plate. 
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To  incorporate  delamination,  the  structure  is  divided  into  regions,  as  seen  in  Fig.  4.1  and 
the  higher  order  displacement  equations  are  applied  to  each  region  (fir  r  =  u,  dl,  d2).  It  must  be 
noted  that  the  refined  higher  order  plate  theory  developed  in  the  previous  chapter  assumed  that 


the  out-of-plane  coordinate,  z,  was  referenced  from  the  midplane  of  the  laminate.  The 
assumption  that  z=0  corresponded  to  the  laminate  midplane  and  that  ±hl2  corresponded  to  the 


surfaces  of  the  plate  allowed  the  simplification  into  the  form  presented  in  Eqs.  (3.32).  Assuming 
that  the  midplane  of  the  original  undelaminated  structure  corresponds  to  z=0,  then  the 
sublaminates  formed  by  the  delamination  are  offset  from  the  original  reference  point.  This 
introduces  eccentricity  in  the  plate  theory  and  Eqs.  (3.32)  need  to  be  modified  to  account  for  this. 
The  displacement  field  of  Eqs.  (3.32)  now  assumes  the  following  form 


(4.1a) 

(4.1b) 

(4.1c) 


where  cr  is  the  coordinate  of  the  local  midplane  of  each  region.  This  introduces  a  change  of 
coordinates  that  permits  the  description  of  the  displacement  that  satisfies  the  conditions  of  zero 


transverse  shear  strains  on  the  upper  and  lower  surfaces  of  the  sublaminate. 


4.1.2.  Continuity  Conditions 

Next,  continuity  conditions  are  imposed  on  the  interfaces  between  the  undelaminated 
region  and  the  sublaminates.  Ideally  it  is  desired  to  match  the  displacement  fields  at  the  interface 
between  the  delaminated  and  undelaminated  regions.  Since  the  displacement  distributions  in  Eqs. 
(4.1)  are  nonlinear,  this  equality  condition  cannot  be  satisfied  exactly.  There  are  a  number  of 
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different  approaches  to  enforce  continuity  between  the  regions.  A  commonly  used  approach  is  to 
satisfy  the  equality  in  an  average  sense.  This  technique  has  the  disadvantage  that  the  slopes  of  the 

out-of-plane  displacement,  —  and  ,  are  generally  not  equal  other  either  side  of  the  interface. 

dx  dy 

In  this  research,  a  more  direct  approach  is  developed  in  which  continuity  conditions  are  imposed 
at  the  interfaces  between  the  undelaminated  region  and  the  sublammates  based  on  the  offset 
between  the  finite  element  nodes  at  the  midplane  of  each  sublaminate.  A  transformation  matrix 
method18  is  used  in  which  the  out-of-plane  displacements  and  slopes  and  the  in-plane 
displacements  are  forced  to  be  equal  for  the  laminate  and  sublaminate  on  each  side  of  the 
interface  resulting  in  the  following  set  of  constraint  equations 

(4 -2a) 
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where  i  represents  each  sublaminate  in  thew  delaminated  region.  A  transformation  matrix  is 
constructed  in  order  to  eliminate  the  degrees  of  freedom  at  the  boundary  of  the  sublammates  from 
the  system  matrices.  This  transformation  condenses  the  system  along  the  boundary  of  the 
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delamination  to  effectively  remove  the  sublaminate  nodes  and  reform  the  system  equations  in 
terms  of  the  nodal  degrees  of  freedom  belonging  to  the  undelaminated  region. 

4.1.3.  Delamination  Crack  Tip  Singularity 

The  formulation  presented  in  the  previous  sections  neglects  the  singularity  created  by  the 
delamination  crack  tip.  As  a  result,  the  strains  and  stresses  in  the  vicinity  of  where  the 
sublaminates  join  the  undelaminated  structure  are  not  accurately  modeled.  Since  the  objective  of 
this  work  is  to  capture  the  effect  of  delamination  on  the  global  response  of  the  structure  this  is  not 
necessarily  of  concern,  provided  that  this  singularity  does  not  create  any  large  changes  in  the 
global  deformation  of  the  structure. 

To  verify  that  the  crack  tip  singularity  has  minimal  effect  on  the  structural  deformation, 
an  idealized  two-dimensional  delaminated  beam  is  modeled.  The  beam,  as  shown  in  Fig.  4.2,  is 
cantilevered  and  has  a  delamination  starting  at  the  free  end.  The  beam  is  assumed  to  be  made  of 
an  isotropic  material  to  allow  convenient  finite  element  modeling.  The  beam  is  modeled  using  a 
traditional  mesh  of  2-D  plane  strain  elements  and,  also,  with  a  2-D  mesh  that  includes  a  finer 
mesh  size  and  quarter-point  crack  tip  elements  in  the  vicinity  of  the  delamination  tip.  The 
quarter-point  crack  tip  elements  are  further  described  in  the  next  chapter.  The  sublaminates  are 
subjected  to  an  opening  force  at  the  free  end,  and  the  resulting  displacement  (end  opening)  is 
computed  for  both  finite  element  meshes.  In  the  following  results,  the  material  is  assumed  to 


Fig.  4.2.  Cantilevered  beam  for  crack  tip  analysis. 
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have  an  elastic  modulus  of  91GPa,  a  Poisson’s  ratio  of  0.3,  and  a  density  of  2000kg/m3.  Figure 

4.3  shows  the  two  finite  element  meshes  in  the  deformed  state. 

Figure  4.4  presents  the  difference  in  results,  between  the  two  models,  for  varying 
delamination  length.  In  this  case,  the  sublaminates  are  each  1mm  thick,  and  the  length  of  the 
undelaminated  section  is  1cm.  In  Fig.  4.4,  the  abscissa  is  the  ratio  of  the  length  of  the 
delamination  to  the  total  beam  thickness.  It  can  be  seen  that  as  the  delamination  becomes  larger, 
the  crack  tip  singularity  has  less  influence  on  the  deformed  shape  of  the  beam. 


Fig.  4.3.  Deformed  shape  of  the  finite  element  meshes  with  and  without  the  crack  tip  singularity. 
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If  the  location  of  the  delamination  within  the  beam  is  changed  to  be  closer  to  the  surface, 
similar  results  are  obtained.  Figure  4.5  shows  the  difference  in  delamination  opening  for  three 
different  delamination  locations.  In  the  figure,  tl/t2,  refers  to  the  ratio  of  thickness  of  the  upper 
sublaminate  to  the  thickness  of  the  lower  sublaminate.  It  is  important  to  note  that  the 
delamination  location  has  virtually  no  effect  on  the  calculated  error.  Since  plate  structures  are 
generally  thinner  than  beams,  neglecting  the  crack  tip  singularity  is  not  expected  to  create  errors 
greater  than  one  percent  for  all  but  the  smallest  of  delaminations.  Small  delaminations  would 
require  a  finer  mesh  of  plate  elements  to  model  the  small  sublaminates,  and  a  more  accurate  field 
description,  through  the  thickness,  may  be  necessary. 


Fig.  4.4.  Error  from  neglect  of  crack  tip  singularity. 
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Fig.  4.5.  Error  from  neglect  of  crack  tip  singularity  for  varying  delamination  locations. 


4.2.  Dynamic  Contact 

The  formulation  for  the  modeling  of  delamination  developed  in  the  previous  section  does 
not  take  into  account  any  possible  contact  between  the  sublaminates.  During  vibration  the 
delaminated  region  opens  up  a  little,  forming  a  gap  between  the  sublaminates  and  then  the 
sublaminates  come  back  together,  impacting  each  other.  Traditional  mode  shapes,  derived  from 
the  eigen  vectors  of  the  plate,  are  not  valid  representations  of  such  behavior  of  a  delaminated 
plate  since  any  opening  between  the  sublaminates  would  correspond  to  a  condition  where  the 
sublaminates  penetrate  each  other  during  the  opposite  half-cycle.  In  order  to  capture  the  true 
response  of  a  plate  during  vibration,  it  is  important  to  model  the  possibility  of  contact  between  the 
sublaminates  in  the  delaminated  region.  Since  little  work  has  been  done  to  assess  the  importance 
of  sublaminate  contact  on  the  accuracy  of  delamination  results,  it  is  included  in  this  model. 

The  contact  between  the  sublaminates  makes  the  response  a  nonlinear  problem  and  time- 
integration  methods  are  generally  used  to  study  the  dynamic  response.  However,  the  sudden 
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impact  of  the  sublaminates  creates  oscillations  in  analytical  response  that  will  destabilize  most 
numerical  methods.  A  discontinuous  time-integration  method  proposed  by  Cho  and  Kim 
(1999a,b)  is  used  in  this  work  to  treat  the  contact  between  the  sublaminates.  This  method  allows 
discontinuity  to  be  incorporated  into  time-integration  without  having  any  other  constraints 
enforced  on  the  nodes  during  contact  and  release,  and  it  also  has  been  shown  to  be  more  stable 
than  other  methods. 

4.2.1.  Contact  Stiffness 

To  represent  contact  between  the  sublaminates,  each  node  of  the  upper  sublaminate  is 
assumed  to  be  connected  to  the  corresponding  node  of  the  lower  sublaminate  by  a  fictitious 
spring.  This  spring  is  assumed  to  have  zero  stiffness  under  tension  and  stiffness  proportional  to 
the  transverse  Young’s  modulus  of  the  plate  when  subjected  to  compression.  This  is  a  typical 
method  often  used  by  other  researchers  modeling  the  buckling  of  composite  laminates  with 
delamination  (Kardomateas  and  Schmueser,  1988;  Pavier  and  Clarke,  1996).  The  difference 
during  dynamic  analysis  is  that  the  sudden  change  in  stiffness  of  the  springs  acts  like  an  impact, 
requiring  careful  modeling  to  correctly  capture.  Since  this  stiffness  matrix  is  not  constant,  the 
dynamic  response  now  becomes  a  nonlinear  problem.  During  the  time  integration  solution  this 
matrix  must  be  recomputed  as  the  contact  surface  between  the  sublaminates  changes  over  time. 

The  contact  stiffness  values  must  then  be  added  to  the  existing  system  equations.  For 
analytical  purposes  the  stiffness  matrix  for  the  contact  springs  is  formed  as  an  independent 
matrix,  kcc,  based  on  the  current  nodal  displacements.  This  stiffness  matrix  is  then  added  to  the 
linear  stiffness  matrix  in  Eq.  (3.76)  before  solving  for  the  nodal  displacements. 

ma  +  cv  +  ku  +  k  ctcu  =  f  (4-3) 

As  noted  before,  the  matrix,  k^,  will  change  over  time  as  the  sublaminates  come  in  contact  with 
each  other  and  then  open  up.  The  bimodularity  of  the  contact  springs  makes  the  contact  problem 
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nonlinear  and  a  time-integration  technique  is  used  to  study  the  response  of  the  laminate  when 
contact  between  the  sublaminates  is  allowed. 

4.2.2.  Discontinuous  Time-Integration  Technique 

The  time-integration  technique  used  in  this  paper  is  based  on  the  discontinuous  time- 
integration  method  of  Cho  and  Kim  (1999a).  The  reason  for  using  a  method  such  as  this,  as 
opposed  to  a  traditional  Newmark  method,  is  that  the  sudden  change  in  stiffness  that  occurs 
during  impact  and  release  of  the  sublaminates  creates  jumps  in  the  field  variables  resulting  in 
numerical  oscillations  in  the  analytical  solution.  The  discontinuous  time-integration  method  uses 
a  generalized  derivative  concept  to  model  the  effect  of  sudden  constraints  that  occur  during 
contact/impact  problems.  The  generalized  derivative  allows  discontinuity  to  be  incorporated  into 
time-integration  without  having  any  other  constraints  enforced  on  the  nodes  during  contact  and 

release. 

The  definition  of  the  generalized  derivative  of  distribution  (Cho  and  Kim,  1999a)  is 
constructed  through  an  integration  by  parts  formula.  With  this  procedure,  the  difficulty  of 
differentiating  a  distribution  is  transferred  to  the  differentiation  of  a  test  function.  Integration  by 
parts  gives  the  generalized  relation  between  the  displacement  vector,  u,  and  the  velocity  vector,  v. 

f"‘  / v  dt  =  -  f"  yT u  dt  +  y Tu["  ^  all  yit)  (4-4) 

**n  n 

where  y  is  the  test  function.  In  weighted  residual  form  this  equation  is  written  as 

0=  f"'yT(v-u)^=  |"*Vrv  +  ^TuJ/-rT«[”+1  for  all  y{t)  (4.5) 

n 

A  similar  equation  can  be  written  for  the  relation  between  the  velocity  and  the  acceleration. 

t'yrzdt  =  -  t'yTvdt  +  yTv\"'  for  all  y(t)  (4-6) 

**n  n  ^ 
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The  time  step  shown  in  Figure  4.6  includes  jump  conditions  for  both  displacements, 
velocities  and  accelerations  at  time  tn,  but  the  test  function  is  assumed  to  be  continuous  to  enforce 
the  effect  of  the  initial  condition. 

Linear  Lagrange  interpolation  functions  are  used  in  the  time  domain  for  approximating 
the  displacement,  u,  velocity,  v,  acceleration,  a,  and  the  test  function,  y.  The  approximation 
vectors  are  written  for  the  interval  t„<  t  <t„+i  as  follows 


u(')=^o('H +^('K 

(4.7a) 

(4.7b) 

a(0  =  ^o('K+^i('K 

(4.7c) 

rit)=Y0{t)rl+y/l{t)r\ 

(4.7d) 

where 

¥o(t)=^ (4.8a) 
(,)=JZ^L_  (4.8b) 

ln+ 1  Ln 

Substituting  Eqs.  (4.7a-d)  into  Eqs.  (4.5)  and  (4.6)  yields  the  approximated  relations  for 
the  time  derivative. 


Fig.  4.6.  Description  of  the  time  domain  with  jumps  in  the  displacement  variable. 


t 


iT.,0  ,  ...  J.1//1,/  v'Tn°  +l//.///.y^Tu),)if|-XTU 


0  =  £ { f 1  ‘ vj  +  «wi ‘ v‘  +  WoK  <+  YiYoY, 

,= 0  " 

0  -  £  {f”‘  (xwiV  +  m  rlV,  +  (w:Iv:  +  twiV  )*}- 


/  V 


86 

(4.9a) 

(4.9b) 


iVj'n  K+YiYtf'n  *n+YiYoYn  V»+Wn 

The  variables  u,  v,  and  a  contain  discontinuities  at  the  initial  time,  t„,  so  the  vectors  from 
the  previous  time  step  are  included. 

■(<,)= ui.,  v(/„)=vi.,  a((,  )=»:_,  r(t,)=rl 

u(<„,)=u;  v(/„1)=v;  a(',.,)=al  r(t,J=r!  (410) 

The  effect  of  the  initial  condition  is  weakly  imposed  via  u‘n_, ,  v„_,  and  a„_, .  Because 
Eqs.  (4.9a, b)  must  hold  for  all  y  'n ,  it  can  be  written  in  matrix  form  as 

®UBtl=®V^+eUB  (4-lla) 

OVn+1  =<DAn+I  +0V„  (4-llb) 


where 


V.„  = 


r« 

I  v‘ 


Anil  = 


n+1  11 


(4.12) 
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(4.14) 


(4.15) 


Inverting  the  relation  yields  the  solution  for  the  time  step 
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n  +  1 
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=  TV„_,  +T0U„ 

(4.16a) 

=  'FA..,  +  V,V. 

(4.16b) 

where 


6 

-^1 

6 

Ai 

.  2 

Ai 

2  J 

(4.17) 


0  I 
0  I 


(4.18) 


The  next  step  is  to  incorporate  the  time  steps  into  the  system  matrices  in  Eq.  (4.3).  It  is 
important  to  note  that  the  system  vectors  u,  v,  and  a  include  not  only  mechanical  displacements, 
but  also  the  electrical  displacements.  The  governing  equations  must  apply  at  both  the  beginning 
and  end  of  each  time  step.  This  yields  following  form 


m  u 

0  m1  1 1  a! 


■  + 


c°  0 
0  c1 


k  +k 
0 


etc 


0 

k1  +k\ 


ff° 

f1 


(4.19) 


or 

MAn+1  +  CVn+1  +(K  +  Kctc)UB+]  =  Fn+,  (4.20) 

Combining  the  Eqs.  (4.16a,b)  and  Eq.(4.20)  gives  the  governing  equation  for  the 
discontinuous  time  step. 

(MT-'r1  +CT"'  +  K  +  Kctc)lJn+1  =Fn+1  -MA.+>  -CV„+1  (4.21) 

where  Vn+i  and  A„+i  are  the  predictors  for  velocity  and  acceleration. 

Vn+1=-T-'JUn  (4.22a) 

An+i  =  -¥'¥  1  JUn  -  T1  JV„  (4.23b) 

These  equations  can  be  solved  using  a  predictor-corrector  solution  algorithm.  Thus,  for 


each  time  step,  the  following  procedure  is  performed 
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(0  Calculate  the  predictors  in  Eqs.  (4.22a, b) 
(h)  Calculate  the  effective  stiffness  matrix 


Keff  =  M'F'1'F'1  +  CT 1  +K  +  Kclc 

(4.24) 

(iii)  Solve  for  current  displacements 

k,/u.„=f„i-ma„,-cv„, 

(4.25) 

(zv)  Correct  the  velocities  and  accelerations 

V„,='K-'U„1+V,., 

(4.26a) 

A„,  =  >r"r,u„l+A„, 

(4.26b) 

This  solution  method  can  be  computationally  expensive,  especially  for  large  system  matrices,  but 
it  effectively  deals  with  the  discontinuity  created  by  contact  between  the  sublammates  and 
prevents  oscillations  in  the  numeric  results.  The  results  shown  by  Cho  and  Kim  (1999a)  illustrate 
how  much  more  effective  this  method  is  than  traditional  Newmark  methods. 


4.3.  Results 

Verification  of  the  proposed  model  is  performed  using  the  experimental  results  obtained 
by  Shen  and  Grady  (1992).  The  test  specimen  examined  is  a  carbon-epoxy  composite 
cantilevered  plate  with  dimensions,  5  in.  length  and  0.5  in.  width.  The  ply  stacking  sequence  is 
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[0°,90o]2s  and  the  ply  thickness  is  0.005  in.,  resulting  in  a  length  to  thickness  ratio  (L/h)  of  125. 
The  material  properties  for  the  carbon-epoxy  composite  are  listed  in  Table  4.1. 

First  the  natural  frequencies  of  the  plate  are  determined  using  the  eigen  values  of  the 
system  and  compared  with  the  experimentally  determined  natural  frequencies  for  the  first 
vibration  mode.  Both  delaminated  and  undelaminated  plates  are  used  in  this  comparison.  For  the 
delaminated  case,  a  centrally  located  (along  the  longitudinal  axis)  through-the-width  delamination 
with  four  different  ply  locations  is  considered.  Delamination  lengths  of  1  in.,  2  in.,  3  in.  and  4  in. 
are  examined  and  compared  with  the  undelaminated  values.  Thus,  in  Table  4.2,  interface  1 
corresponds  to  a  delamination  located  at  the  midplane  and  interface  4  refers  to  the  case  of  surface 
delamination  with  the  uppermost  ply  being  delaminated. 

The  plate  is  modeled  using  a  finite  element  mesh  of  20x4  elements.  The  refined  higher 
order  theory  results  in  seven  degrees-of-freedom  per  node,  or  735  degrees-of-freedom  for  the 
undelaminated  plate.  The  introduction  of  delamination  adds  additional  degrees  for  freedom 
associated  with  the  new  sublaminate.  For  a  2  in.  delamination  the  total  number  of  degrees-of- 
freedom  for  the  system  is  1050. 


Table  4.2. 

Comparison  of  the  natural  frequencies  (Hz)  for  a  delaminated  composite  plate  predicted  by  the 


present  model  with  experimental  data. 


Experimental 

Present  Model 

Undelaminated 

79.83 

82.11 

Interface 

1  in. 

2  in. 

3  in. 

4  in. 

1  -  Exp. 

78.17 

75.38 

66.96 

57.54 

Model 

81.37 

76.73 

67.57 

56.89 

2-  Exp. 

77.79 

75.13 

67.96 

48.37 

Model 

81.42 

77.11 

68.42 

58.04 

3-  Exp. 

80.12 

79.75 

76.96 

72.46 

Model 

81.92 

80.65 

77.44 

71.66 

4-  Exp. 

75.96 

68.92 

62.50 

55.63 

Model 

81.94 

80.83 

78.04 

73.21 
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The  results  from  the  present  model  and  the  experimentally  determined  average  values  are 
listed  in  Table  4.2.  It  is  important  to  note  that  in  the  eigen  analysis,  contact  between  the 
sublaminates  during  vibration  is  not  modeled.  Thus  the  analytically  computed  natural 
frequencies  are  expected  to  over-predict  the  reduction  in  natural  frequency  caused  by 
delamination.  However,  good  overall  correlation  is  observed  between  the  results  from  the  present 
model  and  experiments,  as  shown  in  Table  4.2.  The  only  exception  is  the  case  with  delamination 
located  at  interface  4  (uppermost  lamina).  The  large  drop  in  natural  frequency  that  was  observed 
experimentally  is  most  likely  due  to  buckling  of  the  upper  sublaminate.  To  capture  this 
phenomenon  analytically,  nonlinear  modeling  with  Von  Karman  nonlinear  strains  is  required. 

Next  the  effect  of  contact  between  the  laminates  is  examined.  The  plate  is  subjected  to  a 
one  Newton  tip  load  and  released  from  this  static  position.  The  discontinuous  time  integration 
method  is  used  to  examine  the  effect  of  delamination  on  the  response.  Figure  4.7  shows  the 
analytically  computed  time  histories  of  the  undelaminated  plate  and  the  plate  with  a  2  in. 
delamination  at  interface  2.  The  deflection  plotted  in  Fig.  4.7  is  the  vertical  displacement  at  the 
midpoint  of  the  delaminated  area  on  the  upper  sublaminate.  As  seen  in  Fig.  4.7,  the  response 
curves  for  the  delaminated  plate  with  and  without  contact  are  virtually  superimposed  showing 
that  the  contact  has  little  effect  on  this  low  frequency  vibration  mode.  Figure  4.8  shows  the 
response  of  the  two  sublaminates  for  a  plate  with  a  2  in.  delamination  at  interface  3.  The  opening 
and  closing  of  the  delamination  can  be  clearly  seen,  however,  this  has  little  effect  on  the 
frequency  of  vibration.  These  results  seem  to  indicate  that  the  contact  phenomenon  has  little 
effect  on  the  natural  frequencies  of  low  vibration  modes,  at  least  in  cases  when  the  size  of  the 
delamination  is  not  disproportionately  large.  As  a  result  of  this,  the  natural  frequencies  estimated 
using  the  eigen  values  of  the  delaminated  plate  are  reasonably  accurate,  as  seen  in  Table  4.2.  The 
first  bending  and  twisting  mode  shapes  for  the  undelaminated  plate  are  shown  in  Fig.  4.9,  and 
those  for  the  upper  surface  of  the  plate  with  a  2  in.  delamination  at  interface  2  are  shown  in  Fig. 
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4.10.  The  bending  mode  shows  little  distortion,  but  the  twisting  mode  shows  a  large  amount  of 
twist  concentrated  in  the  location  of  the  delamination. 


Fig.  4.7.  Time  history  for  a  delaminated  plate  under  free  response  from  a  static  load. 


Fig.  4.8.  Displacement  of  upper  and  lower  sublaminates  for  a  delaminated  [0o,90o]2S  plate  under 


free  response. 
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It  is  also  useful  to  examine  the  effects  of  varying  ply  stacking  sequence  on  the  natural 
frequencies  of  the  delaminated  plate.  Figure  4.1 1  shows  the  natural  frequencies  for  the  first  two 
bending  modes  and  the  first  twisting  mode  for  the  plate  with  a  [0o,90o]2S  stacking  sequence  and  a 
2  in.  delamination  located  at  different  interfaces.  It  is  clear  that  delamination  located  towards  the 
midplane  of  the  plate  have  the  largest  effect  on  the  natural  frequency  of  the  first  vibration  mode 
while  delamination  closer  to  the  surface  more  significantly  affects  the  second  and  higher 


Fig.  4.9.  First  bending  and  twisting  modes  for  the  undelaminated  [0°,90°]2S  cantilevered  plate. 


Fig.  4.10.  First  bending  and  twisting  modes  for  the  [0°,90o]2s  cantilevered  plate  with  2  in. 


delamination. 
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frequency  bending  modes.  Also,  larger  impact  of  delamination  is  observed  in  the  twisting  modes. 
Figure  4.12  presents  the  results  for  a  [450,-45°]2S  configuration.  It  can  be  seen  that  the  bending 
modes  of  this  angle  ply  layup  are  less  affected  by  delamination  compared  with  the  crossply  layup. 
This  is  the  result  of  the  increased  anisotropy  of  the  angle  ply  layup  and  the  increased  coupled 
between  the  shear  and  tensile  strains.  However,  the  twisting  frequency  is  reduced  by  over  thirty 
percent.  This  is  an  even  greater  reduction  than  the  twenty-three  percent  drop  observed  in  the 
[0°,90o]2s  laminate. 

So  far  only  balanced  laminates  were  considered.  To  study  the  effect  of  an  unbalance  ply 
layup,  a  laminate  of  [0°,45o]2s  stacking  sequence  is  considered,  and  results  are  shown  in  Fig.  5.13. 
Note  that  although  this  ply  arrangement  causes  significant  coupling  between  the  bending  and 
twisting  modes,  the  reductions  in  natural  frequency  due  to  delamination  are  very  similar  to  those 
observed  in  the  crossply  layup.  Thus,  it  appears  that  the  0°  plies  dominate  the  response  of  the 
cantilevered  composite  plate. 

Next  the  effect  of  laminate  thickness  is  investigated.  The  [0°,90°]2s  laminate  is 
considered,  but  with  an  assumed  ply  thickness  of  0.025  in.,  making  the  laminate  five  times 
thicker  and  changing  the  length  to  thickness  ratio,  L/h  =  25.  In  practice  this  could  be  achieved  by 
stacking  the  plies  in  groups  of  five  to  create  a  [05°,905°]2s  laminate.  The  effects  of  delamination 
on  natural  frequencies  for  this  laminate  are  show  in  Fig.  4.14.  The  natural  frequencies  for  this 
laminate  are  predictably  greater  than  those  of  the  original  laminate  due  to  the  increase  plate 
thickness,  but  the  overall  trend  is  very  similar  to  that  seen  in  Fig.  4.1 1.  The  change  in  thickness 
does  create  moderate  changes  in  the  natural  frequencies  of  the  delaminated  plate,  due  to  the 
significant  transverse  shear  within  the  plate. 
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a)  First  bending  mode 

Fig.  4.14.  Natural 


b)  Second  bending  mode  c)  First  twisting  mode 

frequencies  for  a  cantilevered  [0,9032s  with  L/h  =  25. 


Fig.  4.15.  Displacement  of  the  sublaminates  during  impulse  loading  for  a  plate  with  a  2  in. 


delamination  at  interface  3. 

Finally,  the  response  of  an  adaptive  laminate  in  the  presence  of  delamination  is 
examined.  A  1  in.  by  0.5  in.  piezoelectric  sensor  is  assumed  to  be  bonded  to  the  upper  surface  of 
the  plate.  The  piezo  is  assumed  to  be  made  of  PZT-5H  with  a  thickness  of  0.25  mm  and  is 
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Fig.  4.16.  Comparison  of  PZT  output  voltage  during  impulse  loading  for  a  plate  with  a  2  in. 

delamination  at  interface  3. 

located  at  a  distance  of  0.25  in.  from  the  root  of  the  cantilevered  plate.  The  response  of  the  plate 
is  computed  for  a  5  ps,  IN  impulse  tip  load. 

The  first  case  examined  is  a  [0°)90°]2s  laminate  with  a  2  in.  delamination  at  interface  3. 
Figure  4.15  shows  the  displacement  of  the  upper  and  lower  sublaminates  at  the  center  of  the 
delaminated  zone.  Contact  and  impact  between  the  sublaminates  is  very  visible  due  to  the  high 
frequency  aspects  of  the  vibration.  The  PZT  is  assumed  to  be  open  circuited  and  the  output 
voltage  is  measured.  The  voltage  history  is  shown  in  Fig.  4.16.  It  can  be  seen  that  the 
piezoelectric  voltage  is  definitely  affected  by  the  presence  of  delamination,  but  the  overall  shape 
of  the  output  does  not  change  significantly.  This  implies  that  if  the  PZT  were  used  as  sensor,  the 
output  signal  would  be  similarly  altered  by  the  delamination. 
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Fig.  4.17.  Displacement  of  the  sublaminates  during  impulse  loading  for  a  [0°,90°]2s  laminate  with 

a  2  in.  delamination  at  interface  2. 

The  next  case  analyzed  is  the  same  [0°,90°]2s  laminate,  but  with  the  2  in.  delamination  at 
interface  2.  Figure  4.17  shows  the  displacement  of  the  upper  and  lower  sublaminates  at  the  center 
of  the  delaminated  zone  and  Fig.  4.18  shows  the  voltage  history.  For  this  case,  the  gap  opening 
between  the  sublaminates  is  almost  imperceptible,  but  the  delamination  creates  much  more 
significant  changes  in  the  voltage  output  of  the  sensor.  This  is  consistent  with  the  changes  in 
natural  frequencies  that  were  observed  in  Fig.  4.11.  Little  opening  between  the  sublaminates  is 
also  expected  since  each  sublaminate  is  thick  enough  to  have  appreciable  bending  stiffness. 

The  cases  studied  illustrate  the  difficulty  of  using  low  frequency  response  to  detect 
damage  in  composite  laminates.  Delamination  and  sublaminate  contact  are  more  influential  at 
higher  frequencies,  which  can  be  modeled  using  the  present  approach  with  a  more  refine  finite 
element  mesh.  The  developed  model  can  prove  to  be  a  useful  tool  in  damage  detection  studies. 


PZT  Volta 
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Time  (ms) 


Fig.  4.18.  Comparison  of  PZT  output  voltage  during  impulse  loading  for  a  [0°,90°]2s  laminate 


with  a  2  in.  delamination  at  interface  2. 


5.  Effect  of  Matrix  Cracking 


Unlike  other  studies  that  attempt  to  predict  the  amount  of  damage  based  on  loading 
history,  in  this  study,  it  is  assumed  that  the  damage  is  known.  The  cracking  can  be  in  any 
combination  of  layers  including  the  surface  layers  of  the  laminate.  The  effects  of  matrix  cracking 
must  be  developed  to  reflect  the  structural  model  used.  If  only  the  effects  of  Mode  I  and  Mode  III 
crack  opening  are  considered,  then  the  results  would  then  only  be  appropriate  for  classical  plate 
models  where  transverse  shear  is  not  considered.  For  moderately  thick  plates  modeled  using  a 
first  order  shear  deformation  theory,  the  transverse  shear  plays  a  role  in  plate  deflection,  and 
therefore  Mode  II  crack  behavior  must  be  determined. 

A  number  of  assumptions  must  be  made  in  order  to  model  the  effect  of  matrix  cracking  in 
a  composite  laminate.  The  first  assumption  is  that  the  cracks  can  be  modeled  as  a  statistically 
uniform  array  spaced  at  some  average  distance.  It  is  also  assumed  that  all  of  the  cracks  extend 
vertically  through  the  entire  layer  and  run  parallel  to  the  fibers,  as  shown  in  Figure  5.1.  The  term 
“layer”  is  used  in  this  section  as  opposed  to  “ply.”  The  reason  for  this  is  that  consecutive  plies 
oriented  in  the  same  direction  have  no  means  of  arresting  crack  growth,  which  is  normally  done 
when  the  crack  tip  intersects  fibers  running  in  a  different  direction.  Thus,  multiple  plies  with  the 
same  fiber  orientation  angle  behave  as  a  single,  thick  layer.  As  will  be  described  later,  it  is  also 
assumed  that  cracking  in  one  layer  is  not  influence  by  the  presence  of  cracks  in  other  layers. 
Crack  growth  is  not  considered  in  this  work. 

5.1.  Matrix  Crack  Modeling 

Since  the  equations  of  motion  for  the  composite  plate  developed  earlier  in  this  work  are 
based  on  the  variation  of  the  total  energy,  the  effect  of  matrix  cracking  is  determined  using  the 
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Fig.  5.1.  Assumed  geometry  for  cracked  composite  plate. 


energy  contained  in  crack  opening.  If  a  convenient  method  can  be  found  for  determining  this 
energy  then  the  effect  of  matrix  cracking  can  be  incorporated  as  a  reduction  in  laminate  stiffness. 


n  =  4,5 


(5.1b) 
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where  z  is  the  location  relative  to  the  midplane  of  the  laminate,  and  e°,  k°  and  k2  are 
coefficients  expressed  as  functions  of  w,  v,  w,  \{/x,  and  y/y.  The  total  potential  strain  energy  of 
an  uncracked  plate  can  be  expressed  as  follows 

^ uncracked 

where,  for  convenience,  the  strains  are  separated  into  in-plane  and  out-of-plane  components 

»  =  {4  4  4  <4  <  4  <4  <4  <5-3l> 

4  <  «?r  (5-3b) 

and  the  stiffness  matrices,  Ci  and  C0,  have  been  integrated  through  the  thickness  of  the 
laminate,  as  defined  in  Eqs.  (3.51  and  3.53).  This  notation  allows  convenient  expression 
of  the  stress  strain  relations  in  matrix  form  and  separation  of  the  in-plane  and  out-of¬ 


Vc.s  +  jsXs, 


\dA 


(5.2) 


plane  effects. 

For  a  laminate  with  matrix  cracking  the  total  potential  energy  is  reduced  since  a  portion 
of  the  energy  is  associated  with  the  opening  of  the  cracks  during  deformation.  Thus,  the  total 
potential  energy  is  expressed  as  follows 


U. 


cracked 


-I 


Vc.s  +  isXs, 


dA-AU. 


crack 


(5.4) 


where  the  energy  associated  with  crack  opening,  assuming  that  the  work  done  in  one  cracked 
layer  is  not  affected  by  cracking  in  other  layers  is  defined  as 


Wcrack 


tkpK<  P 


r  Tk 

o  o 


(5.5) 


The  summation  in  this  expression  provides  the  total  energy  for  all  cracked  layers.  The  energy  is 
seen  to  be  proportional  to  the  thickness  of  the  cracked  layer,  tk,  and  the  normalized  crack  density 
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parameter,  //.  The  crack  density  is  defined  as  the  ratio  of  the  thickness  of  the  cracked  layer  to  the 
average  spacing  between  cracks,  ct . 


(5.6) 


The  parameters,  /  and  rak,  represent  the  in-plane  and  out-of-plane  stresses  that  would  be  present 


in  the  uncracked  laminate. 

Computing  crack  opening  energy  at  every  stage  of  the  finite  element  process  would  be  an 
awkward  and  time  consuming  procedure,  so  it  is  desired  to  find  and  an  effective  pair  of  stiffness 

matrices,  Ci  and  C„,  which  have  the  equivalent  energy  as  the  cracked  laminate.  This  then 
satisfies  the  following  expression 


U. 


cracked 


=  1 


— sTCis  +— sJC0s0 

2  2 


— stCjS+— s^C0s0 

2  2  0  0  0 


1  ,Tj 


dA-AU. 


crack 


(5.7) 


The  effect  of  matrix  cracking  can  be  seen  as  a  reduction  in  the  values  of  the  matrix  C, 


yielding  a  new  matrix  C .  It  is  not  sufficient  to  simply  estimate  the  reduction  in  Q  for  the  cracked 

layer  and  apply  this  to  calculate  C ,  because  matrix  crack  opening  causes  dissimilar  reduction  in 
the  extensional  and  bending  stiffnesses.  Therefor,  the  effects  of  matrix  cracking  must  be 
calculated  separately  for  each  element  of  C. 


Fig.  5.2.  Higher  order  stress  distribution  on  the  crack  face. 
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To  determine  the  effects  of  matrix  cracking,  a  procedure  similar  to  that  developed  by 
Gudmundson  and  Zang  (1993)  and  Adolfsson  and  Gudmundson  (1997)  is  used.  But  unlike  their 
work,  which  relies  on  fracture  mechanics  and  the  assumption  of  the  crack  being  in  an  infinite 
homogeneous  medium,  the  current  technique  uses  a  finite  element  approach  based  on  the  actual 
laminate  configuration.  Also,  the  model  developed  by  Adolfsson  and  Gudmundson  (1997)  is 
based  on  the  classical  plate  theory  that  ignores  transverse  shear  stresses.  Since  a  higher  order 
model  is  used  in  this  paper,  additional  parameters  need  to  be  determined. 

Using  the  principle  of  superposition,  the  effective  stresses  in  the  cracked  laminate  are 
equal  to  the  summation  of  the  stresses  corresponding  to  c2,  cr4,  and  a6  in  the  local  coordinate 
system,  and  to  modes  I,  II,  and  m  crack  loading,  as  seen  in  Figure  5.2.  For  the  kth  layer,  this 
loading  can  easily  be  separated  into  components  as  follows 

T  =  t0+t]tj  +  t2t]2 +riT]3  for  modes  I  and  HI  (5.8) 

=  *oo  +  To^  +  To2r?2  for  mode  E  (5-9) 


where 


(5.10) 


where  zk  is  the  distance  of  the  layer  midplane  from  the  laminate  midplane.  Thus  the  stresses  on 
the  crack  face,  t  and  xQ,  are  cubic  and  quadratic  functions  respectively.  In  matrix  form  the 
stresses  are  expressed  as 


(5.11a) 


t 
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(5.11b) 

These  components  are  defined  based  on  the  strains  in  the  uncracked  laminate  as  follows 


(5.12a) 

T,  =|tc,(^  +  3z,V) 

(5.12b) 

=-^”Tcl  zkK2 

(5.12c) 

^4Tc'^ 

(5.1 2d) 

r..=T.c  .(f.°+z,V) 

(5.13a) 

!■„  =',T„c„  zhK; 

(5.13b) 

(5.13c) 

where  T  and  T0  are  the  coordinate  transformation  matrices  and  c,  and  c0  are  the  stress-strain 

relations  for  in-plane  and  out-of-plane  components. 

The  next  step  is  to  relate  the  stresses  to  the  global  strains.  The  global  stress  vectors 

(5 . 1 1  a,b)  are  rearranged  to  form  the  following  relations 

I3  Zk\  3  Zk^  3 

0  |i,  !<,*,% 

0  0  ^t\zk  I3 

0  0  ^-I3 

8  3 


Tc,  0  0  ' 

0  Tc,  0  s  (5.14a) 

0  0  Tc, 
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^2  Zk\l 

®  ^kZiP  2 


0 


I, 


T0c0  0 
.  0  T0cOJ 


(5.14b) 


where  Ij  is  an  identity  matrix  of  rank  j.  Substituting  Eqs.  (5.14a,b)  into  Eq.  (5.5)  yields  the 
following 


At/„,  =  E 


k  =  \ 
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’TC< 

0 

1 

O 

Tc, 

0 

0  " 

sT 
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Tc, 

0 

pkk 
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Tc, 

0 

1  A  k 

T*  P 

1 

0 

0 

Tcj 

0 

0 

Tc,_ 

+  s 


T0c0  0 
0  T0c0 


T0c0  0 
.  0  T0cOJ 


dA 


(5.15) 


In  the  local  coordinate  system  the  matrices,  and  Plf  ,  are  diagonal  matrices  relating 
the  applied  stress  to  crack  opening  energy. 


P1*  = 


kk 

01 

kk 


nkk 

HQ! 


fi* 

PS  Pn 


PI o 

nkk  nkk  nkk 

PlQ 


pr  = 


>25  P° o“ 
K  P°? 


i  j 


(5.16a) 


(5.16b) 


The  component  matrices  are  defined  for  each  geometric  variation  of  the  stress  as  follows 


Pt  = 


0  0 
i  k 


0 


0  flc.) 

0  0  p* 


(3) 


(5.17a) 


Pof  = 


0*  0 

(2) 

0  0 


(5.17b) 


where  the  nonzero  components,  Pi  {V) .  P\  (2)  anc*  P\  (3)  >  are  stress  to  ener6y  relations  for  the 
specific  crack  opening  modes,  I,  II  and  HI.  The  other  components  are  all  zero  since  the  matrices 
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are 


defined  based  on  the  local  coordinate  system.  It  can  also  be  shown  (Gudmundson  and  Zang, 
1993)  that  for  an  elastic  material  the  component  matrices  are  symmetric  and  that 

P?=Pf  and  y 6of=(3of  (5-18) 

This  results  in  a  system  of  equations  with  fifteen  unknown  parameters.  Combining  Eqs. 
(5.7)  and  (5.15)  gives  the  reduced  stiffness  matrices. 


N 

-z 

Tc, 

0 

0  " 

Tc, 

0 

0  ' 

c,  =c, 

tkpk 

0 

Tc, 

0 

0“ 

0 

Tc, 

0 

(5.19a) 

k= 1 
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Co 
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-t\ 

tk  pk 

1  \ 

•  c 

o 

0 

T.c._ 

Pf 

T.c. 

0 

0 

TX.. 

1 

(5.19b) 

Once  the  damage  parameters  are  known  substitution  into  Eqs.  (5.19a,b)  gives  the  stiffness 
matrices  for  the  laminate  containing  matrix  cracking.  It  should  be  noted  that  this  formulation 
based  on  crack  tip  energy  does  not  take  into  account  crack  growth,  thus  it  is  not  useful  for 
predicting  further  damage  evolution.  The  next  step  will  be  to  develop  an  effective  method  for 
calculation  the  damage  parameters  for  an  arbitrary  laminate. 


5.1.2.  Finite  Element  Analysis  of  a  Representative  Crack 

These  parameters  could  be  estimated  using  closed  form  solutions  if  assumptions  like 
those  of  Adolfsson  and  Gudmundson  (1997)  are  made.  Adolfsson  and  Gudmundson  (1997)  used 
fracture  mechanics  solutions  based  on  stress  intensity  factors  to  relate  the  load  to  the  work  done. 
Because  of  the  difficulty  in  obtaining  fracture  mechanics  solutions,  it  was  necessary  to  use  the 
solutions  for  cracks  in  infinite  homogeneous  media.  The  accuracy  of  this  assumption  varies 
dramatically  based  on  laminate  thickness  and  material  properties. 

To  develop  a  more  reliable  approach,  a  separate  finite  element  analysis  is  used  to  relate 
the  individual  loads  to  changes  in  work.  The  separate  finite  element  modeling  of  the  crack 


Material  coordinate 
axes  for  cracked 
layer 

►  2 


Fig.  5.3.  Representative  matrix  crack  used  to  compute  damage  parameters. 


behavior  in  each  layer  with  cracks  causes  definite  increases  in  the  computational  effort. 
However,  the  developed  procedure  offers  advantages  that  can  offset  this  increase  in  CPU  time. 
First,  the  (3  matrix  can  be  derived  for  a  number  of  different  crack  densities,  then  a  simple  curve 
fitting  process  can  be  used  to  express  p  as  a  function  of  crack  density.  Second,  this  function  can 
then  be  stored  and  reused  if  the  same  laminate  configuration  is  encountered  again.  Therefor, 
though  computational  effort  is  required  the  first  time  a  new  material  or  laminate  configuration  is 
used,  subsequent  analyses  will  require  very  minimal  calculations. 

The  method  developed  in  this  work  takes  a  single  crack  and  meshes  it  with  a  2-D  finite 
element  grid.  The  method  then  enforces  the  condition  of  zero  global  strain.  Stresses  based  on  the 
components  of  global  strain  are  then  applied  sequentially  to  the  crack  face.  The  work  done  by  the 
stress  is  then  calculated  for  each  case.  Using  these  values  of  work,  the  damage  parameters  for 
Eqs.  (5.17a, b)  are  then  computed.  The  details  of  this  procedure  will  next  be  explained. 
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Fiber  direction 
in  the  cracked 
layer 


Fig.  5.4.  Geometry  of  perpendicular  plane  used  for  the  finite  element  analysis. 


Fig.  5.5.  Finite  element  mesh  for  representative  crack. 

Based  on  the  assumption  of  a  uniform  distribution  of  matrix  cracks  a  single  cracked 
segment  of  the  laminate  is  considered,  as  shown  in  Figure  5.3.  The  cross  section  perpendicular  to 
the  length  of  the  crack,  shown  in  Figure  5.4,  is  meshed  with  a  two-dimensional  finite  element 
grid.  Because  the  matrix  crack  runs  parallel  to  the  fibers,  the  cross  section  is  oriented 
perpendicular  to  the  fibers  and  is  in  the  material  2-3  plane.  The  use  of  this  cross  section  orients 
the  coordinate  system  so  that  no  coordinate  transformation  is  necessary  during  these 
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computations.  The  resulting  mesh  is  shown  in  Fig.  5.5  with  the  deformation  induced  by  the 
applied  stresses. 

The  finite  elements  are  eight  noded  serendipity  quadrilateral  elements  with  three  degrees 
of  freedom  per  node,  as  shown  in  Fig.  5.6a.  Two  degrees  of  freedom  represent  the  usual 
displacements  within  the  plane  of  the  mesh,  but  a  third  degree  of  freedom  is  added  to  represent 
out-of-plane  displacement  caused  by  pure  shear  (o,2  in  the  material  coordinate  system  or  Mode 
HI  crack  opening).  These  displacements  are  referred  to  as  uh  u2  and  u3,  corresponding  to 
displacements  along  the  material  1-,  2-  and  3-  axes.  The  nodes  on  the  sides  of  the  elements  are 
located  at  the  midpoint  for  all  elements  except  those  at  the  crack  tip.  The  elements  adjacent  to  the 
crack  tip  are  quarter  point  crack  tip  elements,  as  seen  in  Fig.  5.6b,  with  nodes  at  the  quarter  point 
from  the  tip  in  order  to  model  the  singularity  that  exists  at  the  crack  tip.  This  method  ensures  that 

accurate  results  are  obtained  using  a  reasonable  number  of  nodes. 

To  obtain  the  three-dimensional  stress  tensor  from  the  two-dimensional  mesh,  the 
assumption  is  made  that  the  displacements  do  not  vary  with  respect  to  the  coordinate 
perpendicular  to  the  plane  of  the  mesh.  Since  the  2-D  mesh  is  aligned  with  the  2-3  material  axes, 
the  following  assumption  is  made. 

=  =  =  0  (5.20) 

etc,  9x,  fix, 


The  strain  is  now  approximated  as  follows 


This  strain  vector  is  then  used  to  compute  the  elemental  stiffness  using  the  finite  element  method. 

The  element  stiffnesses  are  computed  for  each  element  and  added  to  the  global  system 
matrix.  Since  the  work  must  be  computed  for  crack  opening  under  zero  global  strain,  boundary 
conditions  are  enforced  on  the  nodes  on  the  end  of  the  finite  element  mesh.  Also,  because  the 
crack  section  is  from  the  middle  of  a  distribution  of  cracks,  the  ends  must  represent  a  symmetric 
boundary  with  other  crack  sections.  Based  on  these  two  conditions  the  displacement  for  all  nodes 
on  the  end  are  fixed  in  the  material  1-  and  2-  directions  and  the  displacement  of  the  midpoint  of 
the  end  is  fixed  in  the  3-  direction.  The  system  stiffness  matrix  is  symmetric  positive  definite,  so 
a  Cholevsky  decomposition  can  be  used  to  solve  the  system. 


Fig.  5.6.  Serendipity  (a)  and  crack  tip  (b)  finite  elements  used  to  analyze  representative  matrix 


crack. 
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There  are  fifteen  unknown  parameters,  so  fifteen  work  computations  must  be  made  to 
gain  enough  information  to  solve  for  every  parameter.  Crack  opening  occurs  only  for  those 
strains  which  create  a  stress  on  the  crack  surface.  Thus  the  only  strains  that  need  to  be  considered 
are  combinations  of  s\,  k\,  k\,  £°,  k°6,  k26,  e°4,  k\  .  The  stress  applied  to  the  crack  face  is 
then  calculated  using  the  following 


V 

0  ,1 

<r2 

>  =  <v 

£°  +z(lC2  +Z2Kj  ) 

£°  +  z(k°  +  Z2K26  )j 

(5.22a) 


(5.22b) 


Fifteen  combinations  of  the  strains  s\ ,  /c°,  k\,  £°,  fC°,  k\ ,  £4 ,  k4  are  used  to  create  fifteen 

stress  distributions  on  the  crack  face.  The  deformation  associated  with  each  combination  is  then 
computed  using  back-substitution  on  the  decomposed  stiffness  matrix.  The  work  is  then 
computed  by  integrating  the  product  of  the  stress  and  displacement  over  the  crack  face. 

Using  Eq.  (5.15)  the  work  for  the  following  fifteen  combinations  of  strain  can  be 


computed 
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(5,23b) 

(5.23c) 
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k°2  and  k\ 

=>  W6~~  (C22  )2  |^00(1)  (K2  )  +  Pum  {K2  )  +  2fiom)K2  K2  ] 

•  (5.23f) 
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|  (5.23p) 

Using  the  work  computed  using  the  finite  element  model  for  each  case,  the  matrix 


cracking  damage  parameters  can  be  computed. 
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(5.24c) 

(5.24d) 
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This  procedure  is  conducted  for  each  layer  containing  matrix  cracking.  Once  the  damage 
parameters  have  been  computed,  they  can  be  used  give  the  reduced  laminate  stiffness  matrix 
when  analyzing  the  response  of  the  structure  containing  the  cracking. 


5.2.  Bimodularity  Effects 

Modeling  of  a  composite  plate  with  matrix  cracking  and  delamination  will  be 
accomplished  using  the  following  steps.  First,  reduced  stiffness  values  will  be  determined  for 
layers  containing  matrix  cracks  using  the  separate  finite  element  model  described  above.  Each 
cracked  layer  will  be  analyzed  individually  and  given  its  own  reduced  stiffness  values.  The  plate 
will  then  be  modeled  using  the  higher  order  theory  and  the  finite  element  method.  Once  the  finite 
element  stiffness  matrix  for  the  plate  is  formed,  continuity  conditions  will  be  enforced  between 
the  regions  as  described  in  the  previous  chapter.  However,  a  problem  that  is  encountered  is  that 
Mode  I  crack  opening  does  not  occur  when  the  crack  is  under  compression.  Thus,  the  composite 
laminate  has  a  different  stiffness  depending  upon  the  state  of  stress.  This  phenomenon  is  known 
as  bimodularity. 

For  static  problems,  bimodularity  can  be  addressed  by  using  an  iterative  technique.  The 
neutral  axis,  the  location  in  the  z  direction  where  the  stress  is  zero,  is  first  assumed  and 
compressive  stiffness  properties  are  applied  to  the  layers  under  compression  and  tensile  stiffness 
properties  are  applied  to  the  layers  under  tension.  The  displacements  are  then  determined,  and  a 
new  neutral  axis  is  calculated.  The  stiffness  matrix  is  recalculated  using  the  new  value  and  the 
problem  is  solved  again.  The  iteration  continues  until  the  neutral  axis  converges  to  a  final  value. 
For  the  case  of  matrix  cracking,  few  iterations  are  required  since  the  difference  in  overall 
laminate  stiffness  for  varying  locations  of  the  neutral  axis  is  relatively  small,  compared  to  the 
bimodularity  that  results  in  such  materials  as  Kevlar-rubber  composites. 
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Fig.  5.7.  Natural  frequency  for  mode  2  inflection  point  location  for  [0,902]2s  glass-epoxy 

laminate. 

When  solving  for  plate  vibration,  the  nonlinearity  associated  with  the  bimodularity  of 
matrix  cracking  means  that  superposition  of  vibration  modes  is  not  applicable.  Natural  frequency 
values  can  no  longer  be  calculated  using  standard  eigenvalue  analysis  since  the  plate  stiffness  is 
dependent  upon  the  deflection  shape,  or  more  specifically,  whether  the  plate  is  in  compression  or 
tension.  The  natural  frequencies  can  be  solved  for  individually  if  the  mode  shape  was  known  and 
used  to  calculate  the  correct  stiffness  matrix.  Fortunately,  the  bimodularity  introduced  by  matrix 
cracking  is  relatively  small  and  results  in  only  very  small  changes  in  the  mode  shape  from  the 
uncracked  case  to  the  case  with  matrix  cracks.  To  demonstrate  this,  the  technique  of  Wu,  Jing 
and  Chin  (1989)  is  used  to  solve  a  transfer  matrix  model  of  a  cantilever  beam  with  matrix  cracks. 
The  natural  frequency  is  calculated  assuming  different  locations  for  the  inflection  points,  and  the 
local  maximum  or  minimum  represents  the  actual  natural  frequency.  The  results  are  shown  for  a 
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[0,902]2s  glass-epoxy  laminate  with  saturation  cracking  in  all  of  the  90°  layers  (Pk  1.0).  This 
model  was  chosen  since  matrix  cracking  has  a  large  effect  in  glass-epoxy  and  laminates  with  a 
large  number  of  90°  plies.  Figure  5.7  shows  the  results  corresponding  to  the  second  mode  of 
vibration.  The  symbol  “0”  represents  the  inflection  point  for  the  uncracked  laminate  and  it  is 
apparent  that  the  inflection  point  has  shifted  only  a  very  small  distance.  Also,  the  error 
introduced  in  the  natural  frequency,  if  the  calculations  were  made  using  the  uncracked  mode 
shape,  is  on  the  order  of  0.01%.  This  is  so  small  because  near  the  inflection  point,  the  bending 
stresses  are  always  small  and  small  changes  in  the  stiffness  have  little  effect.  The  third  mode, 
shown  in  Fig.  5.8,  demonstrates  similar  behavior.  However,  in  this  case,  there  are  two  inflection 
points,  resulting  in  a  surface  plot.  Again,  the  shift  in  inflection  point  locations  is  small  and  has 
little  effect  on  the  calculated  natural  frequency. 


Fig.  5.8.  Natural  frequency  for  mode  3  inflection  point  location  for  [0,902]2s  glass-epoxy 


laminate. 
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The  results  from  this  model  show  that  good  estimates  of  the  natural  frequency  can  be 
obtained  by  using  the  uncracked  mode  shapes  in  the  calculation  of  the  stiffness  matrix.  The 
accuracy  is  further  improved  by  using  an  iterative  process  during  the  natural  frequency 
calculation  where  the  uncracked  mode  shape  is  used  as  an  initial  starting  point  in  the  calculation 
of  the  stiffness  matrix.  The  natural  frequency  and  mode  shape  are  determined  next  using  eigen 
value  analysis  and  the  new  mode  shape  is  used  to  recalculate  the  stiffness  matrix.  This  process  is 
continued  until  convergence  of  the  inflection  points  and  neutral  axes  are  achieved. 

During  plate  vibration,  bimodularity  can  create  further  complexities  if  varying  amounts 
of  cracking  exist  in  each  of  the  different  layers.  In  such  case,  the  structure  will  have  different 
stiffness  depending  upon  whether  the  plate  is  being  bent  up  or  down.  Thus,  the  period  for  each 
half-cycle  during  vibration  at  a  natural  frequency  is  different,  and  the  overall  natural  frequency 

(to)  is  determined  by 

2 

CO  =  - 17  (5.25) 

+  CO  2 

where  0)i  and  co2  are  the  natural  frequencies  for  each  half-cycle  (Bert,  et  al.,  1989).  In  other 
words,  to  calculate  the  overall  natural  frequency  (co)  the  plate  is  analyzed  first  for  the  upward 
deflection  case  (to,),  and  then  again  for  the  downward  deflection  case  (co2). 

5.3.  Results 

In  this  section  results  are  presented  using  the  developed  model  for  transverse  matrix 
cracking.  The  reduction  in  laminate  stiffness  caused  by  matrix  cracking  is  first  examined  and  the 
results  from  the  model  are  compared  with  experimental  data  and  results  from  other  researchers. 
Next,  matrix  cracking  is  included  in  the  modeling  vibration  in  laminate  structures,  and  its  effect 


on  the  dynamic  response  is  examined. 
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Table  5.1. 


Material  properties  used  for  matrix  cracking  comparison. 

Property 

Graphite-epoxy 

Glass-epoxy 

E,  (GPa) 

144.8 

41.7 

E2  (GPa) 

9.6 

13.0 

Vl2 

0.31 

0.30 

v23 

0.461 

0.42 

G12  (GPa) 

4.8 

3.40 

G23  (GPa) 

3.3 

4.58 

5.3.1.  Stiffness  Reduction 

The  first  step  is  to  verify  that  the  developed  model  accurately  predicts  the  loss  of 
laminate  stiffness.  This  is  done  by  modeling  sample  composite  plates  and  comparing  the  results 
with  published  experimental  results.  Results  show  that  this  technique  accurately  predicts  the 
reduction  in  stiffness  for  a  variety  of  materials  and  laminate  configurations. 

Shown  in  Figs.  5.9  and  5.10  are  the  reductions  in  extensional  and  bending  stiffness  for 
cracked  laminates  of  [0,902]s  graphite-epoxy  and  [0,903]s  glass-epoxy,  respectively.  The  results 
obtained  from  the  present  model  are  compared  with  the  experimental  results  from  Highsmith  and 
Reifsnider  (1982)  for  the  glass-epoxy  laminate.  For  the  graphite-epoxy  laminate  the  experimental 
work  of  Groves  (1986)  is  used.  The  material  properties  for  the  two  composite  laminates  are  listed 
in  Table  5.1.  Though  the  elastic  moduli  of  the  two  composite  materials  are  very  different,  the 
present  model  shows  good  correlation  in  both  cases.  This  is  in  contrast  to  other  models  whose 
accuracy  depends  upon  the  type  of  composite  laminate  examined. 

Due  to  a  lack  of  sufficient  experimental  data  on  the  effects  of  bending  in  the  literature, 
ANSYS  finite  element  models  are  used  for  validation.  Finite  element  modeling  of  a  cracked 
structure  is  difficult  to  accomplish  effectively  and  efficiently,  thus  motivating  the  current 
research.  Individual  matrix  cracks  are  spaced  apart  by  a  distance  on  the  order  of  one  or  two  ply 
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thicknesses,  thus  a  single  composite  specimen  will  likely  have  hundreds  of  matrix  cracks.  The 
desire  for  accurate  results  motivates  the  modeler  to  use  a  fairly  large  number  of  elements  in  the 
vicinity  of  each  crack  and  particularly  near  the  tips.  The  geometry  of  the  crack  and  the  variety  of 
loading  possibilities  motivates  the  use  of  a  3D  analysis.  Such  a  model  would  result  in  a  very 
large  number  of  degrees  of  freedom.  Thus  for  verification  two  slightly  simplified  NASTRAN 
and  ANSYS  models  are  used.  The  first  is  a  2D  plane-strain  NASTRAN  model  of  the  cross 
section  of  a  composite  laminate  containing  a  small  number  of  matrix  cracks.  The  spacing  of  the 
cracks  is  varied  slightly  to  reflect  a  realistic  distribution.  The  plate  model  is  then  loaded  and  the 
average  strains  are  compared  to  the  results  from  an  uncracked  model  of  the  same  composite  cross 
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Crack  Density  (t/a) 

Fig.  5.10.  Stiffness  loss  of  glass-epoxy  [0,903]s  laminate. 


section.  The  second  model  is  a  3D  ANSYS  mesh  of  a  single  representative  crack.  The  rational 
behind  this  model  is  a  comparison  versus  the  plane  strain  assumption.  The  size  of  this  model  is 
considerably  lager  in  spite  of  modeling  a  single  crack,  and  enforcement  of  the  boundary 
conditions  involved  significant  effort.  The  reduction  in  bending  stiffness  is  significantly  smaller 
compared  to  the  reduction  in  extensional  stiffness,  due  to  the  fact  that  matrix  crack  opening  is 
much  smaller  during  bending  than  it  is  during  extension.  Due  to  this  reason,  models  that  are 
designed  for  predicting  extensional  stiffness  should  not  be  used  to  determine  bending  stiffness. 

Figure  5.11  presents  a  comparison  of  the  present  model  with  the  experimental  results 
from  Highsmith  and  Reifsnider  (1982),  as  well  as  with  the  other  available  techniques.  These 
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techniques  include  the  work  by  Allen  and  Lee  (1991),  Hashin  (1986)  and  Adolfsson  and 
Gudmundson  (1997).  The  present  model  correlates  well  with  the  experimental  results  and 
asymptotically  approaches  the  ply  discount  method.  Unlike  Hashin’ s  method,  which  can  only  be 
used  to  study  a  few  specialized  cross-ply  laminates,  the  present  approach  is  applicable  to  any 
laminate  configuration. 

The  results  shown  thus  far  have  all  examined  the  extensional  stiffness  of  cross-ply 
laminates.  Matrix  cracking  also  affects  the  shear  stiffness  of  the  laminate.  This  affect  can  be 
seen  by  not  only  looking  at  the  shear  stiffness  of  cross-ply  laminates,  but  also  by  examining  the 
extensional  stiffness  of  angle  ply  laminates.  Figure  5.12  shows  the  reduction  in  stiffness  for  a 
[0,±45]s  glass-epoxy  laminate.  In  this  case,  the  crack  occurs  only  in  the  angle  ply  layers  and  both 


Fig.  5.11.  Normalized  extensional  stiffness  of  glass-epoxy  [0,903]s  laminate. 
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the  +45°  and  -45°  layers  are  assumed  to  have  equal  crack  densities.  Again  the  current  model  has 
good  agreement  with  the  experimental  data  (Highsmith,  et  al.,  1981),  even  in  the  presence  of 
increased  in-plane  shear  forces  due  to  the  angle  plies.  Figure  5.13  presents  the  reduction  in  the 
shear  modulus  for  a  [0,902]s  graphite-epoxy  cross-ply  laminate  with  cracking  in  the  inner  layer  as 
compared  to  the  experimental  data  of  Tsai  and  Daniel  (1992).  This  is  the  only  experimental  data 
available  in  the  literature  on  the  affect  on  matrix  cracking  on  shear  stiffness.  Again  the  developed 
model  shows  good  correlation  with  the  experimental  results.  There  is  a  more  noticeable 
difference  between  the  model  and  experiments  than  seen  in  previous  results.  This  is  most  likely 
due  to  the  difficulty  in  accurately  measuring  shear  the  modulus  of  composite  laminates  and  it 
should  be  noted  that  the  experimental  results  do  not  seem  to  converge  to  zero  reduction  for  zero 
crack  density. 


Fig.  5.12.  Extensional  stiffness  loss  of  [0 ,±45],  glass-epoxy  laminate. 
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Fig.  5.13.  Reduction  in  shear  modulus  of  [0,902]s  graphite-epoxy  laminate. 

5.3.2 .  Plate  Vibration  with  Matrix  Cracking 

Figure  5.14  shows  the  effect  of  matrix  crack  density  on  natural  frequency  for  the  same 
[0,902] 2s  glass-epoxy  cantilever  plate  discussed  above.  The  first  case  is  for  cracking  only  in  the 
upper  90°  layer  and  the  second  case  is  for  cracking  in  all  90°  layers.  The  asymmetric  amount  of 
cracking  in  case  1  results  in  differing  frequencies  for  each  half-cycle  (coi  and  co2  in  Fig.  5.14). 
This  effect  is  most  apparent  for  mode  1,  where  all  of  the  cracks  are  closed  during  the  upward 
deflection,  coj.  Since  the  cracks  close,  their  only  effect  is  on  transverse  shear,  and  they  have  little 
effect  on  natural  frequency.  During  the  downward  deflection,  0)2,  all  of  the  cracks  open  causing  a 
significant  loss  in  stiffness,  thus  reducing  the  natural  frequency.  For  case  2,  since  the  cracking  is 
symmetric,  the  upward  and  downward  deflections  have  the  same  behavior,  and  the  half-cycles  do 
not  have  different  frequencies.  For  the  higher  modes,  the  upward  and  downward  deflections 
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create  both  areas  with  closed  cracks  and  areas  with  open  cracks,  thus  decreasing  the  difference  in 
the  two  half-cycle  frequencies. 

Next,  results  are  presented  for  cases  where  matrix  cracks  are  accompanied  by 
delamination.  Figure  5.15  shows  an  example  of  how  matrix  cracking  can  have  significant  effect 
on  structural  behavior  that  is  often  not  considered.  Shown  in  this  figure  is  the  effect  of 
delamination  in  a  [0,90]2s  graphite-epoxy  cantilever  plate,  as  measured  by  Shen  and  Grady 
(1992).  The  example  shown  is  for  a  through  the  width  delamination  between  the  first  90°  layer 
and  the  second  0°  layer.  Results  are  shown  for  the  case  of  no  matrix  cracks  and  the  case  of 
cracking  in  all  90°  layers  (p=1.0).  The  present  theory  uses  a  finite  element  mesh  with  40 
elements  and  504  degrees  of  freedom.  It  can  be  seen  that  the  effect  of  matrix  cracking  is 
comparable  in  magnitude  to  the  effect  of  the  delamination.  Also,  delamination  area  increases  the 
effect  of  matrix  cracking,  as  seen  by  a  larger  reduction  in  natural  frequency  at  high  delamination 


Fig.  5.15.  Natural  frequency  of  a  delaminated  graphite-epoxy  [0,902]2s  cantilever  plate. 
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lengths.  This  is  due  to  the  fact  that  the  delamination  creates  a  new  surface  within  the  laminate 
and  cracks  originally  internal  to  the  laminate  are  now  on  the  surface,  which  dramatically 
increases  their  effect  on  stiffiiess.  This  example  shows  that  matrix  cracking  and  delamination 
need  to  be  addressed  simultaneously  when  studying  the  effect  of  damage  on  composite  laminates. 


6.  Adaptive  Structure  Optimization 


Having  modeled  the  adaptive  structural  system,  the  next  consideration  is  how  to  optimize 
that  system  for  maximum  performance.  Optimization  of  adaptive  structures  is  a  popular  area  of 
study,  but  only  a  limited  amount  of  work  has  been  done  on  simultaneously  optimizing  both 
structural  and  electrical  aspects  of  a  system.  A  common  application  of  the  electrical  interaction 
with  the  structural  deformation  is  in  the  design  of  passive  electrical  damping  systems.  Work  has 
been  done  on  defining  optimal  electrical  parameters  for  the  damping  circuit,  but,  in  general,  very 
simple  structural  models  were  used  that  do  not  take  into  account  the  complex  state  of  strain  that 
may  exist  in  the  piezoelectric  patches.  In  addition,  the  placement  of  the  piezoelectric  patch  has 
not  been  considered  concurrently  with  the  design  of  the  electrical  system. 

In  this  chapter  it  is  demonstrated  how  multidisciplinary  optimization  (MDO)  techniques 
can  be  utilized  to  optimize  both  structural  and  electrical  aspects  of  an  adaptive  structural  system. 
To  simultaneously  optimize  multiple  performance  requirements,  the  Kreisselmeier-Steinhauser 
(K-S)  function  approach  (Wren,  1989;  Sethi  and  Striz,  1997)  is  used.  This  technique  combines 
all  the  objective  functions  and  the  constraints  to  form  a  single  unconstrained  composite  function 
to  be  minimized.  Then  an  unconstrained  solver  is  used  to  locate  the  minimum  of  the  composite 
function.  Optimization  of  an  integrated  structural  and  electrical  system  involves  both  discrete  and 
continuous  variables.  The  optimization  method  used  in  this  work  is  very  similar  to  the  hybrid 
method  of  Seeley,  Chattopadhyay  and  Brei  (1996),  with  the  major  difference  being  that  the 
discrete  search  and  the  gradient  based  search  are  more  independent  due  to  the  nature  of  this 
particular  problem.  Since,  the  present  model  assumes  that  only  the  electrical  components  to  be 
continuous  parameters  the  gradient  based  search  essentially  amounts  to  finding  the  optimal  set  of 
electrical  components  for  a  give  structural  configuration.  Composite  plates  of  varying  stacking 
sequence,  boundary  conditions  and  surface  bonded  actuators  are  considered. 
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6. 1.  Problem  Formulation  /Description 

The  objective  of  the  optimization  procedure  is  to  minimize  the  peak  response  of  a 
composite  plate  with  piezoelectric  patches  connected  to  a  passive  electrical  damping  circuit.  The 
optimization  procedure  must  determine  the  optimal  combination  of  plate  geometry  and  electrical 
components  that  will  result  in  a  tuned  damping  circuit  with  the  minimum  response  for  a  given  set 
of  modes.  The  passive  electrical  damping  circuit  is  assumed  to  be  composed  of  linear  inductors, 
resistors  and  capacitors  with  values  that  are  determined  during  the  optimization  process.  Note 
that  the  optimization  procedure  only  determines  the  optimal  values  for  the  components  and  not 
the  configuration  of  the  damping  circuit.  Other  nonlinear  electrical  components  could  be  used 
and  optimized,  but  a  corresponding  increase  in  computational  effort  would  be  expected.  The 
optimization  procedure  must  also  take  into  account  location  and  orientation  of  the  piezoelectric 
patches  on  the  plate  and  determine  the  location  that  provides  for  maximum  damping. 

The  mathematical  optimization  problem  is  stated  as  follows 

Min  f(<j>)  (6.1) 

where 

■  xc„ 

YCh 

ANGh 

*=  R, 

Li 
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h  =  l...Na 

i  =  l...Nr 
j  =  l...Nj 
k  =  1...N. 
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subject  to  the  constraints 
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where  f(<t>)  is  the  objective  function  representing  the  peak  response  of  the  system  for  a  set  of 
vibrational  modes,  §  is  the  vector  of  design  variables,  Na  is  the  number  of  actuators  and  Nr,  Nj, 
and  Nc  are  the  numbers  of  resistors,  inductors  and  capacitors  to  be  optimized.  The  design 
variables  include  the  x-  coordinate  (XC),  the  y-  coordinate  (YC)  and  the  orientation  angle  (ANG) 
for  each  actuator  and  the  resistance  (R),  inductance  (L)  and  the  capacitance  values  for  each 
electrical  component.  It  is  also  possible  to  include  ply  orientation  angles  as  design  variables. 
Geometric  constraints  are  imposed  so  that  the  electrical  components  maintain  positive  values  in 
order  to  represent  physical  hardware. 

6.2.  Optimization  Technique 

The  model  for  piezoelectric  system  developed  in  previous  chapters  is  used  to  model  the 
structural  system.  For  this  analysis  the  structural  design  variables  (XC,  YC,  ANG)  are  all  chosen 
to  be  discrete  variables,  while  the  electrical  design  variables  (R,  L,  C)  are  all  continuous.  The  x- 
and  y-  coordinates  of  the  piezoelectric  patch  are  chosen  as  discrete  variables  based  on  the  finite 
element  model  of  the  plate  structure.  The  plate  is  meshed  with  a  uniform  grid  of  elements  and  the 
piezoelectric  patches  are  required  to  have  locations  that  align  the  PZT  with  the  mesh.  This 
procedure  is  used  to  avoid  remeshing  of  the  plate  during  every  optimization  iteration.  The  small 
numerical  variations  that  result  from  changing  the  mesh  in  a  finite  element  analysis  can  introduce 
errors  that  may  lead  to  sub-optimal  designs.  Thus,  it  is  more  efficient  to  use  a  uniform  mesh  and 
restrict  the  PZTs  to  discrete  locations.  The  analysis  can  start  with  a  relatively  coarse  mesh  to 
determine  the  general  coordinates  of  the  optimum  location,  and  then  if  further  accuracy  is  desired, 
the  process  can  be  repeated  using  a  more  refined  mesh  and  the  previous  optimum  solution  as  a 
starting  point.  The  orientation  angles  are  included  to  allow  modeling  of  piezoelectric  patches  that 
are  not  isotropic  within  the  plane  of  the  patch  and  exert  an  actuation  force  in  a  preferential 
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direction.  These  include  inter-digitated  electrode  (IDE)  piezoelectric  patches  (Rodgers,  Bent  and 
Hagood,  1996)  which  are  very  powerful  since  they  make  use  of  d33  actuation. 

The  optimization  formulation  is  based  on  the  Kreisselmeier-Steinhauser  (K-S)  function 
approach.  In  this  technique  the  multiple  objective  functions  are  first  transformed  into  reduced 
objective  functions,  which  can  be  expressed  as  follows 

f.  _  (p.  (<j))/Fi0  )-l-gmax  <0  i  =  1 . . .  number  of  objective  functions  (6.4) 

where  Fi0  represents  the  original  value  of  the  ith  objective  function,  F;  is  its  value  based  on  the 
current  design  variables  and  gmax  is  the  largest  value  of  the  original  constraint  vector.  These 
normalized  objective  functions  are  now  analogous  to  constraints  and  are  combined  into  a  single 
vector,  fm(<|>);  m=l,2,...,M,  with  M  being  the  sum  of  the  number  of  constraints  and  the  number  of 
objective  functions.  The  constraints  and  objective  functions  are  then  combined  into  a  single 
composite  objective  function,  FksW?  defined  as 

FKsW  =  f«x+1taIeP(f"(*H'“)  <«> 

P  m=l 

where  f^  is  the  largest  constraint  corresponding  to  the  new  constraint  vector  fm(<|>).  The 
parameter  p  acts  as  a  draw-down  factor  controlling  the  distance  from  the  surface  of  the  K-S 

envelope  to  the  surface  of  the  maximum  constraint  function. 

The  optimization  algorithm  consists  of  a  hybrid  scheme  that  uses  simulated  annealing  to 
optimize  the  discrete  variables  and  a  continuous  search  procedure  to  optimize  the  continuous 
variables.  Because  the  K-S  function  approach  reduces  the  problem  to  an  unconstrained  one,  any 
unconstrained  gradient  based  technique  can  be  used  for  the  continuous  variables.  A  modified 
Newton  method  is  used  for  the  unconstrained.  A  finite  difference  scheme  is  used  to  calculate  the 
gradient  and  the  Hessian  matrix,  necessary  for  determining  the  search  direction.  The 
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computational  cost  associated  with  the  second-order  search  is  mostly  overcome  by  convergence 
in  significantly  fewer  iterations  compared  to  other  first  and  zero  order  search  techniques. 

The  optimization  algorithm  begins  with  the  simulated  annealing  procedure  and  a  user 
defined  initial  point.  The  simulated  annealing  procedure  controls  the  values  of  the  discrete 
variables  only,  and  attempts  to  minimize  the  objective  function  using  a  probabilistic  approach. 
Whenever  the  value  of  the  objective  function  is  requested  by  the  simulated  annealing  algorithm, 
the  current  values  of  the  discrete  variables  are  used  to  calculate  the  open  circuit  eigen  values  and 
eigen  vectors  of  the  plate.  Then  an  unconstrained  search  procedure  is  invoked  to  minimize  the 
objective  function  at  that  point  by  varying  the  continuous  design  variables.  Since  the  continuous 
variables  are  the  electrical  components,  the  procedure  can  be  described  as  determining  the 
optimum  passive  damping  circuit  for  a  particular  set  of  PZT  locations  and  orientations.  The 
objective  function  value  returned  to  the  simulated  annealing  search  is  the  minimum  achievable  K- 
S  function  value  determined  by  the  unconstrained  search. 

When  the  value  of  the  objective  function  is  required  by  the  unconstrained  search,  the 
current  values  for  the  continuous  variables  are  used  to  construct  the  linear  electrical  matrices. 
These  are  then  combined  with  the  reduced  structural  matrices  based  on  the  eigen  values  and 
vectors.  These  equations  are  then  solved  to  obtain  the  system  response  for  any  disturbance 
frequency.  Next,  the  frequency  response  curves  are  calculated  and  the  maximum  response  peak 
is  computed  at  each  mode  to  be  minimized.  These  parameters  are  then  used  to  calculate  the  K-S 
function  value. 


Fig.  6.1 .  Configuration  for  the  cantilevered  plate  prior  to  optimization. 


132 


V 


6.3.  Results  for  optimal  passive  damping  systems 

The  developed  optimization  procedure  is  demonstrated  first  through  the  integrated  design 
of  a  cantilevered  plate  with  a  single  actuator  connected  to  a  passive  shunt  circuit.  The 
cantilevered  plate  has  been  studied  extensively  in  the  literature  related  to  passive  damping  and 
therefor  is  a  good  example  for  benchmarking  the  developed  optimization  algorithm.  This  is  due 
to  the  fact  that  for  all  bending  modes,  the  root  is  subjected  to  maximum  strain,  thus  making  it  the 
optimal  location  for  a  piezoelectric  actuator  to  generate  the  maximum  possibly  polarization 
energy  and  thereby  controlling  vibration. 

The  plate  is  assumed  to  be  made  of  3.17mm  thick  aluminum  and  the  detailed  dimensions 
are  shown  in  Figure  6.1.  A  single  piezoelectric  actuator  is  placed  on  one  side  to  induce  vibration 
and  a  single  actuator  is  placed  on  the  opposite  side  and  is  connected  to  the  shunt  circuit.  The 
plate  is  modeled  with  a  19x3  element  finite  element  mesh  and  the  system  is  reduced  using  the 
first  twelve  modes.  Optimization  is  performed  to  minimize  the  vibratory  response  associated 
with  the  first  mode.  The  design  variables  include  location  of  the  piezoelectric  patch  and  the 
electrical  parameters  governing  the  inductor  and  the  resistor  in  the  passive  shunt.  In  the  initial 
design,  the  patch  is  assumed  to  be  at  the  center  of  the  plate.  Results  are  presented  for  two 
configurations,  one  with  the  inductor  and  resistor  in  series  and  the  other  in  parallel. 

The  frequency  response  curves  before  and  after  optimization  are  shown  in  Figure  6.2.  It 
can  be  seen  that  the  optimization  algorithm  is  able  to  reduce  the  peak  response  of  the  first  mode 
by  over  an  order  of  magnitude  in  both  cases.  The  optimum  location  for  the  piezoelectric  patch  is 
determined  by  the  algorithm  to  be  at  the  root  of  the  plate,  as  was  expected  for  a  cantilevered  plate. 
The  optimum  values  calculated  for  the  inductor  and  resistor  in  series  are  422.032  henries  and 
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Fig.  6.2.  Frequency  response  curves  for  a  cantilever  plate  with  and  without  passive  damping 


circuits. 

14406  ohms,  respectively.  For  the  case  with  the  inductor  and  resistor  in  parallel,  the 

corresponding  values  are  417.927  henries  and  489630  ohms. 

Next  a  carbon  fiber-epoxy  composite  plate  simply  supported  on  all  four  sides  is 
considered.  This  example  provides  a  more  interesting  design  challenge  due  to  the  complexity  of 
the  mode  shapes,  which  are  associated  with  different  locations  and  orientations  of  maximum 
strain.  In  such  cases,  it  is  possible  that  locations  and  orientations  that  are  optimal  for  a  particular 
mode  may  fall  on  the  inflection  lines  where  strain  is  minimal  for  the  other  modes.  In  such  a  case 
the  piezoelectric  actuator  would  not  generate  any  polarization  energy  during  vibration  that  excites 
that  particular  mode,  making  passive  damping  of  that  mode  impossible.  The  piezoelectric 
actuators  used  in  this  example  are  Active-Fiber  Composite  (AFC)  (Rodgers,  Bent  and  Hagood, 
1996)  actuators  to  allow  orientation  angle  of  the  actuators  to  be  included  as  design  variables.  The 
plate  dimensions  are  32cm  by  32cm  with  1 .6mm  thickness.  The  initial  lay-up  for  the  plate  is 
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eight  laminae  in  [0,90]2s  configuration.  Vibration  is  induced  by  a  2cmx2cm  piezoelectric  actuator 
located  on  the  backside  of  the  plate.  The  plate  is  modeled  with  a  16x16  finite  element  mesh  with 
a  single  4cmx4cm  AFC  actuator  located  in  the  center  of  the  plate  at  zero  orientation  angle  as  an 
initial  design. 

First,  the  AFC  patch  is  modeled  with  a  parallel  shunt  circuit  containing  one  inductor  and 
one  resistor  in  parallel.  The  system  is  then  optimized  to  reduce  the  response  peak  for  each  of  the 
first  four  vibrational  modes.  The  frequency  response  curves  associated  with  the  initial  and  the 
optimum  configurations  are  presented  in  Figures  6.3-6.  In  each  case  the  targeted  mode  is  reduced 
by  more  than  a  factor  of  ten,  showing  the  effectiveness  of  the  developed  optimization  algorithm 
in  damping  each  of  the  individual  modes.  The  optimized  location  of  the  piezoelectric  patch,  for 
the  first  mode,  is  the  comer  of  the  plate,  as  illustrated  in  Figure  6.7.  The  optimized  orientation 
angle  for  the  AFC  patch  is  45°  and  the  inductor  and  resistor  assumes  values  of  110.645  henries 
and  558307  ohms,  respectively.  For  the  second  and  third  modes,  the  optimal  location  is  predicted 
to  be  a  quarter  of  the  way  inwards  from  the  center  of  one  side.  These  locations  coincide  with  the 
center  of  the  maximum  out-of-plane  deflection  as  shown  in  Figures  6.8  and  6.9.  For  the  second 
mode  the  optimal  actuator  angle  is  90°,  and  the  inductor  and  resistor  have  values  of  20.8077 
henries  and  177846  ohms,  respectively.  For  the  third  mode  the  optimal  actuator  angle  is  0°,  and 
the  inductor  and  resistor  have  values  of  12.2968  henries  and  138118  ohms,  respectively.  The 
optimal  location  for  the  forth  mode  predicted  by  the  optimization  algorithm  is  shown  in  Figure 
6.10.  The  optimal  angle  for  the  actuator  is  60°,  and  the  inductor  and  resistor  have  values  of 
7.2863  henries  and  192467  ohms,  respectively.  This  location  is  not  intuitively  obvious  to  be  the 
optimal  location  for  the  forth  mode.  However,  convergence  to  this  solution  was  obtained 
irrespective  of  the  initial  design.  This  particular  case  further  illustrates  the  necessity  of  using 
formal  optimization  techniques  in  the  design  of  structures  with  complex  mode  shapes. 
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Fig.  6.5.  Frequency  response  curve  for  the  optimal  passive  damping  design  using  mode  3. 
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Fig.  6.6.  Frequency  response  curve  for  the  optimal  passive  damping  design  using  mode  4 
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Fig.  6.7.  Optimized  actuator  location,  (a),  for  passive  damping  of  the  first  vibration  mode,  (b). 
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Fig.  6.8.  Optimized  actuator  location,  (a),  for  passive  damping  of  the  second  vibration 


mode,  (b). 
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Fig.  6.9.  Optimized  actuator  location,  (a),  for  passive  damping  of  the  third  vibration  mode,  (b). 
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Fig.  6. 10.  Optimized  actuator  location,  (a),  for  passive  damping  of  the  fourth  vibration 


mode,  (b). 
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The  optimum  orientation  angles  of  the  AFC,  for  each  case  studied,  are  different.  This  is 
because,  in  each  case,  the  piezoelectric  fibers  of  the  AFC  actuators  align  with  the  direction  of 
maximum  strain.  The  results  obtained  demonstrate  the  benefit  of  including  both  structural  and 
electrical  parameters  simultaneously  in  the  optimization  procedure.  In  the  numerical  examples 
shown,  optimal  designs  are  obtained  by  not  only  varying  the  parameters  of  the  passive  shunt 
circuit,  but  by  also  tailoring  the  mechanical  stiffiiess  of  the  piezoelectric  patch.  The  location  of 
the  piezoelectric  patch  not  only  affects  the  passive  damping  capability,  it  also  changes  the 
dynamic  response  of  the  system  causing  natural  frequencies  to  shift.  This  may  or  may  not  be 
desirable  for  a  particular  problem.  In  such  cases  the  natural  frequencies  can  be  included  as 

additional  objective  functions  using  the  K-S  function  approach. 

Next  the  same  plate  is  optimized  with  two  2cmx4cm  AFC  actuators.  Each  actuator  is 
connected  to  a  single  parallel  shunt  circuit.  In  this  example  the  objective  is  the  minimization  of 
the  first  three  vibrational  modes.  This  method  is  very  different  from  the  multi-mode  damping 
circuit  proposed  by  Wu  (1998)  which  uses  three  different  circuits  connected  to  a  single 
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Fig.  6.11.  Position  of  the  actuators  on  the  simply  supported  plate  before  and  after 
optimization  for  two  actuators  and  three  modes. 


140 


t 


1.00E-05  ; 


£  1.00E-06 


c 

<D 

E 

01 

o 

JS 

CL 

< n 


g  1.00E-07 


1.00E-08 

0 


100  150  200  250 


300  350  400 


Forcing  Frequency  (Hz) 

Fig.  6.12.  Frequency  response  curves  for  a  simply  supported  plate  both  undamped  and 

with  two  actuators. 


piezoelectric  actuator  in  order  to  damp  the  three  modes.  The  method  developed  by  Wu  would  be 
less  effective  in  this  case  since  the  optimal  locations  for  each  mode  are  different  and  the  use  of  a 
single  actuator  would  require  compromises  in  the  placement  on  the  plate.  The  optimized  results 
of  the  circuit  are  shown  in  Figures  6.11  and  6.12.  The  piezoelectric  patches  move  to  locations 
corresponding  to  5cmxl2cm  and  15cmx22cm,  with  orientations  of  0°  and  75°  respectively,  as 
seen  in  Figure  6.11.  The  inductor  and  resistor  assume  values  of  24.2490  henries  and  620624 
ohms  for  the  first  actuator  and  41.5787  henries  and  531028  ohms  for  the  second  actuator,  after 
optimization.  The  resulting  frequency  response  curve  is  presented  in  Figure  6.12.  It  can  be  seen 
that  this  damping  circuit  affects  only  the  second  and  third  vibration  modes.  This  is  because 
circuits  with  a  single  inductor  and  a  single  resistor  are  effective  over  a  narrow  range  of 
frequencies.  Thus,  in  this  case,  two  actuators  with  single  circuits  cannot  damp  out  more  than  two 
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of  the  first  three  modes.  The  optimal  values  of  the  design  variables  also  roughly  correspond  to  the 
locations  and  angles  predicted  for  the  second  and  third  modes  presented  in  the  single  actuator 
case. 

In  the  final  case  studied,  the  effect  of  composite  tailoring  is  investigated  by  including  the 
ply  orientations  of  the  plate  as  additional  design  variables.  The  same  plate  with  two  2cmx4cm 
AFC  actuators  is  considered.  The  composite  laminae  are  allowed  to  be  oriented  at  even  multiples 
of  15°,  but  symmetry  is  enforced  to  reduce  the  number  of  design  variables.  Thus,  only  four 
additional  design  variables  are  introduced.  The  optimized  results  are  shown  in  Figures  6.13  and 
6.14.  The  piezoelectric  patches  move  to  locations  corresponding  to  15cmxl0cm  and 
21cmxl2cm,  with  orientations  of  75°  and  15°  respectively,  as  seen  in  Figure  6.13.  The  inductor 
and  resistor  have  values  of  38.828  henries  and  559984  ohms  for  the  first  actuator  and  211.850 
henries  and  230981  ohms  for  the  second  actuator.  The  final  ply  lay-up  for  the  composite 
laminate  is  [-15°,  75°,  15°,  -15°]s.  The  resulting  frequency  response  curve  is  shown  in  Figure 
6.14  along  with  the  response  curve  for  the  initial  cross-ply  design.  It  can  be  seen  that  the 
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Fig.  6.13.  Position  of  the  actuators  on  the  simply  supported  plate  before  and  after 
optimization  for  two  actuators  and  three  modes,  including  optimization  of  ply  orientation. 
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Fig.  6.14.  Frequency  response  curves  for  the  undamped  initial  design  of  a  simply  supported  plate 


and  the  optimal  design  including  two  actuators  and  optimized  ply  orientations. 


inclusion  of  the  ply  angles  are  effective  in  reducing  the  overall  response  across  the  frequency 
spectrum  thus  aiding  the  passive  damping  actuators  which  are  effective  at  specific  frequencies. 

The  numerical  examples  presented  demonstrate  the  potential  of  using  MDO  techniques 
for  designing  integrated  adaptive  structural  systems.  For  passive  damping  circuits,  optimization 
determines  best  possible  values  of  both  structural  and  electrical  components  and  can  be  used  to 
design  systems  for  damping  of  multiple  modes.  But  greater  value  lies  in  the  potential  of 
designing  systems  with  synergistic  characteristics,  where  careful  combinations  of  actuators  and 
circuits  provide  significant  damping  with  little  or  no  power  consumption.  Also,  this  technique 
can  model  the  behavior  of  real  electrical  components  as  opposed  to  the  ideal  components 
considered  thus  far.  Factors,  such  as  the  internal  resistance  of  inductors,  can  be  included  in  the 
model  without  having  to  reformulate  the  optimization  algorithm. 


* 


7.  Concluding  Remarks 

A  smart  structural  model  has  been  developed  to  analytically  determine  the  response  of 
arbitrary  structures  with  piezoelectric  materials  and  attached  electrical  circuitry.  The  developed 
model  includes  the  effects  of  composite  damage  in  the  form  of  delamination  and  matrix  cracking. 
The  mathematical  model  uses  a  coupled  piezoelectric-mechanical  theory  that  accurately  captures 
both  electrical  and  mechanical  characteristics  of  adaptive  structures.  Parametric  studies  were 
performed  to  assess  the  difference  in  results  predicted  between  the  refined  higher  order  laminate 
theory  and  the  classical  plate  theory,  as  well  as  between  the  coupled  and  uncoupled  piezoelectric 
models.  The  effect  of  matrix  cracking  is  included  as  a  reduction  in  laminate  stiffness  calculated 
using  a  separate  finite  element  procedure.  Delaminated  plates  are  modeled  as  a  collection  of 
sublaminates  and  contact  between  the  sublaminates  is  modeled  using  a  discontinuous  time 
integration  method.  A  robust  multi-objective  optimization  procedure  was  developed  to  optimally 
determine  both  the  placement  of  piezoelectric  actuators  and  parameters  describing  associated 
electrical  components  in  a  passively  damped  structure.  The  following  conclusions  were  made 
from  the  present  study 

1 .  The  developed  finite  element  model  gives  accurate  results  for  the  response  of  adaptive 
composite  laminates  with  piezoelectric  devices. 

2.  The  model  is  able  to  model  the  interaction  between  a  structure  and  attached  electrical 
circuitry,  and  good  correlation  is  observed  with  experimental  data  for  passive  electrical 
damping  circuits. 

3.  For  thin  plates,  the  differences  in  mechanical  displacements  predicted  by  the  higher  order  and 
the  classical  plate  theories  are  significant  only  when  the  time  step  is  sufficiently  small  to 
capture  very  high  frequency  modes.  The  differences  in  sensor  output  are  more  significant  in 


all  cases. 
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4.  For  thicker  plates,  significant  differences  exist  in  both  mechanical  displacement's  and  sensor 
output  predicted  by  the  higher  order  and  the  classical  plate  theories  during  all  analyses. 

5 .  The  assumption  of  constant  electric  field  over  the  entire  PZT  area  in  the  uncoupled  theory 
leads  to  inaccuracies  in  modeling  the  effect  of  high  frequency  vibrations  that  create  both 
areas  of  local  compression  and  local  tension  within  the  PZT  patch. 

6.  The  natural  frequencies  for  the  low  frequency  modes  are,  in  general,  accurately  estimated 
using  eigen  values  without  contact  effects.  However,  contact  between  the  sublaminates 
becomes  more  important  in  the  presence  of  higher  frequency  vibrations. 

7.  Delamination  embedded  closer  to  the  midplane  of  the  laminate  shows  greater  affect  on  the 
natural  frequencies,  but  is  less  affected  by  contact  between  the  sublammates. 

8.  The  crack  tip  singularity  of  a  delamination  does  not  significantly  affect  the  global  response  of 

the  structure. 

9.  The  method  used  to  model  matrix  cracking  is  able  to  accurately  predict  the  loss  of  stiffness 
for  a  variety  of  laminate  materials  and  ply  layups,  based  on  good  correlation  with 
experimental  data  for  the  change  in  extensional  and  shear  stiffness. 

10.  Matrix  cracking  introduces  nonlinearity  into  a  structure  due  to  crack  closure  under 
compression.  This  effect  can  be  modeled  using  techniques  for  bimodular  materials. 

1 1 .  The  developed  method  is  capable  of  optimizing  passive  damping  parameters  in  various 
electrical  configurations,  as  demonstrated  by  optimization  of  both  parallel  and  series 
configurations. 

12.  The  optimized  locations  of  the  piezoelectric  actuators  correspond  to  the  locations  of 
maximum  strain  for  the  mode  being  optimized,  and  the  fibers  of  Active  Fiber  Composites  are 
oriented  so  that  the  fibers  are  parallel  to  the  principle  strain. 
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